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Abstract 
Introduction. When milling complex-profile surfaces of parts, the selection of tool trajectories and orientations affect 


the roughness parameters. However, in the studies devoted to the formation of trajectories, recommendations to provide 
the quality of microgeometry of surfaces were not taken into account. Moreover, when writing programs for CNC 
equipment in CAM systems, the limitations of cutting modes were determined exclusively using a geometric approach. 
It did not take into account the influence of the orientation angles of the sphero-cylindrical tool relative to the normal 
plane on the quality of surface treatment, namely on roughness. The work was aimed at the creation of the methodology 
for selecting the limiting values of the orientation angles of a sphero-cylindrical tool to optimize the process of 
machining spatially complex surfaces. The tasks included achieving the minimum values of the amplitude roughness 
parameter Rz and determining the effectiveness of various machining paths. 

Materials and Methods. Methods of correlation and regression analysis were used, the results were compared and 
generalized. The least-squares method was applied to estimate the parameters of the regression equation. The DMU 50 
ecoline processing center was used for the experimental studies. Roughness was measured on a Surfcam 1800 D 
profilometer. The material of the samples was steel 12X18N10T. The material of the tool was hard alloy 1620 Sandvik 
with PVD coating (physical vapor deposition, the closest domestic analogue is T15K6). 

Results. It has been shown in detail how roughness parameters Rz depend on the angle of inclination and the diameter 
of the tool. Twenty examples were summarized in a table. Natural regression coefficients were calculated using linear 
and hyperbolic models. It was found that the diameter of the tool had a greater effect on the formation of roughness 
parameter Rz than the angle of inclination. For a detailed description of the influence features, the coefficients of 
multiple, partial, paired correlation and multiple determination were compared. The limitations associated with the 
angles of inclination of the tool when processing complex surfaces were determined. A scheme for calculating the angle 
of the normal was visualized, which included the selected step along the axis to determine the lengths of the segments 
of the broken curve. The profilograms of surfaces obtained with different shaping trajectories were given in the form of 
drawings. This allowed us to conclude that milling from top to bottom is unsuitable when the tool is tilted 5°-35°. A 
map has been compiled by which it is possible to judge the roughness, knowing the type of milling and the inclination 
angle (from 5° to 80°). The dependence of the roughness parameter on the processing speed and the use of coolant was 
represented graphically. The calculated parameters for determining the optimal angle of inclination of the tool were 
tabulated. Their analysis proved the adequacy of the proposed method of preparing control information. 

Discussion and Conclusion. The presented technique made it possible to determine the optimal values of the 
orientation angles of the sphero-cylindrical tool. taking into account the cutting speed and the minimum possible 


amplitude roughness parameter Rz. The pattern of feeding /2 = 0.4 mm/tooth for surface areas with a total angle of 
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5—50° was considered. In this case, processing along trajectories in the passing, opposite and bottom-top directions, 
provided roughness in the range of 3—6 wm according to parameter Rz. The top-down toolpath is not recommended for 


use in final operations due to the significant height of parameter Rz. 


Keywords: amplitude roughness parameter, orientation of a sphero-cylindrical tool, milling of complex-profile 


surfaces, spatially complex surfaces 
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AHHOTalna 

Beedenue. (pu ppesepoBanuu ci0x%KHbIX MOBepXHOCTeH AeTaei BbIOOp TpaeKTopHii HW OpHeHTaljMu MHCTpyMeHtTa 
BJIMAIOT Ha MapaMeTpbl WepoxoBaTOcTH. OWHakO B HCCJICOBaHHAX, MOCBALICHHbIX (POPMHpOBaHHO TpaeKTOpHii, He 
YUHTBIBAIOTCA PCEKOMCH AHH, MO3BOIAIOWIMe OOeCIeYHTh KaYeCTBO MHKpOreoMeTpHH MoBepxHoctTel. K Tomy 2%«Ke mpv 
HallucaHHu TporpaMM Ayia oOopygoBanua c UITY Bp CAM-cuctremax (oT aura. computer-aided manufacturing — 
AaBTOMATH3HPOBaHHOe MpOW3BOACTBO) OFpaHHYeHHA PexXKUMOB Pe3aHHA OMPeCMAIOTCA MCKJIIOUNTEIbHO C TMOMOLIbIO 
reoMeTpHyeckoro Toszxofa. OH He YYNTHIBACT BIMAHHE yIIOB OpHeHTaluu ciPepouMIHHApHyeckoro HHCTpyMeHTa 
OTHOCHTEJIbBHO IWIOCKOCTH HOpMasIM Ha KayecTBO OOpaboTKH HoBepxHocTeli, a HMeHHO Ha wWepoxoBaTocTs. Lemp 
paOoTbI — co3qaHve MeTOAMKM 0 BbIOOpy MpexeIbHbIX 3HaYeHHH yrIOB OpHeHTalMu ccbepouMIMHApHyeckoro 
WMHCTPyMeHTa JIA ONTHMM3allMH Wpolecca MexaHHyecKkol OOpaOoTKH MpOCcTpaHCTBeEHHO-CJIOXKHbIX MOBepXHOcTel. 
3aqauu: JOCTWKeHHe MMHMMAJIbHBIX 3HaYeCHHM aMIWVIMTy{HOrO MapaMeTpa wWepoxoBatoctu Rz u onpeyesenne 
3PbeKTHBHOCTH pa3IHUHbIX TpaeckTopuli OOpaooTKu. 

Mamepuansi u memoooi. VUicnonb30Baiucb MeTObI KOPpeAWMOHHOTO MW perpeccHOHHOTO aHasIv3a, pe3YyJIbTAThI 
cpaBHuBaIMch HM oOobmamMch. Ja OeCHKH MapaMeTpoB ypaBHeHuA perpeccHu MpHMeHAICA MeTO HaMMeHbIIHXx 
KBagpaTosB. Jit 9KCIepHMeHTaJIbHbIX HCCIeqOBAHHU 3afevicTBOBaNM OOpadaTEIBarolMi weHTp DMU 50 ecoline. 
IWepoxosatoctTs u3MepssIH Ha Ipodpusometpe Surfcom 1800 D. Matepuan o6pa3n0B — ctamb 12X18H10T. Matepnan 
MHCTpyMeHTa — TBepybIit criaB 1620 Sandvik c PVD-noxppituem (oT aur. physical vapor deposition — dbu3n4eckoe 
ocaxkeHHe MapoB MeTasWIOB. ONWKaMIIM OTeUeCTBeHHBIM aHaslor — T15K6). 

Pe3yiemamoi uccnedosanua. J\etaibHo 0Ka3aHO, Kak WlapaMeTpbI WepoxoBaTocTH Rz 3aBHCAT OT yrula HaKJIOHa U 
yjMaMetpa WHcTpyMeHta. J[BalljaTb MpHMepoB pefcTaBeHbl B Bue TaOmuuUbI. EcTecTBeHHbIe KOIPMUIMeCHTHI 
perpeccuu paccuuTaHbI 0 JIMHeMHOM u rMnepOomMyecKoH MoyemaM. YcTaHOBJIeHO, YTO AMaMeTp MHCTpyMeHTa 
OolbWIe BIMAeT Ha (OPMHMpOBaHHe MapamMeTpa WepoxoBaToctn Rz, yem yrou HaksoHa. J[a WeTabHOrO omMcaHHA 
ocoOeHHOcTeli BIIMAHHA CpaBHUBaIMch KOIPPHUMCHTbI MHOXKECTBCHHOM, 4aCTHOH, WapHoli KoppesayuH 
MHOXKECTBEHHOM WeTepMuHaHH. OnpeyeseHbl orpaHHyeHHa, CBA3aHHbIe C yrylaMM HaKJIOHa MHCTpyMeHTa TIpv 
oOpaboTke cNO%KHEIX MoBepxHoctei. BusyanM3upoBaHa cxeMa JIA pacueTa yrula HOpMasIM, KOTOpad BKJIIOUaAeT 


BbIOpaHHbIt War WO OCH AIA OMpeAeNeHHA IM OTPe3KOB JOMaHol KpuBOH. JjaHbI B BUC PHCyHKOB MpodusorpaMMbI 
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MOBEPXHOCTeH, MOJYYeHHbIe MPH pa3sIM4HbIX TpaeKTOpHAx PopMoobpa30BaHHA. ITO MO3BONMIO CAelaTb BbIBO O 
HeIIpHrogHOcTH :pesepoBaHuaA cBepxy BHH3 Mp HakIoHe HHCTpyMeHta 5°—-35°. CoctapmeHa KapTa, 10 KOTOpoH 
MOXKHO CYAHTb O WepOxoBaTOCTH, 3Had BU (pexepoBaHHaA M yrou HakoHa (oT 5° yo 80°). Ppaduuecku moKa3ana 
3aBHCHMOCTb MapaMeTpa WepOxoBaTOCTH OT CKOpocTH OOpaboTKH HM MpHMeHeHHA OxIaKAalollel %XUAKOCTH. CBeeHbI 
B TaOMUy pacueTHble MapaMeTpbl JIA ONpeyeeHHA ONTHMaIbHOrO yruia HakIOHa MHCTpyMeHTa. Ux aHaM3 WoKa3zas 
ajleKBaTHOCTb IIpeAIOKCHHOLO MeTOAa NOATOTOBKH yiipaBAOWen HHPopMallHn. 

O@6cyarcdenue u 3akmo4uenue. I[peyctaBneHHad MeTOAUKa MO3BONMIa ONpPeeIMTb OMTHMAIbHbIeC 3HAYCHHA YTJIOB 
opHeHTalyHu cibepouMIMHypMyeckoro WHCTpyMeHTa Cc y4eTOM CKOpOCTH pe3aHHA HM OCTHXKeHHA MUHMMAJIbHO 
BO3MOXHOrO aMIVIHTyIHOrO MapaMetpa wepoxosatoctu Rz. Paccmotpena cutyauua nogaun fz=0.4 mM/3syO wna 
y4acTKOB HOBepXHOCTH C CyMMapHbIM yrsiom 5°-50°. B sTom cyyyae oOpaboTKa 0 TpaeKTOpHAM B TIOILyTHOM, 
BCTpeYHOM HallpaBJIeHHH MU CHH3y BBepx OOecreuHa WepoxoBaTOCcTb B jMala30He 3-6 MKM TIO Mapametpy Rz. 
TpaekTopua JBYwKeHHA CBepXy BHH3 He PeKOMeHJOBaHa K TIPHMeCHCHHIO Ha OKOHUYATECIIBHBIX OlepallHAx 13-3a 


3Ha4UTeIbHOHM BBICOTHI TlapaMeTpa Rz. 


KaroyeBbie CJ10Ba: aMIWIMTYIHbI MapaMeTp WepoxoBaTOcTH, OpHeHtallua ciepouMIMHApHyeckoro HHCTpyMeHTa, 


dbpezsepoBaHue CJIOPKHbIX MOBepxHOcTeH, TIPOCTpaHCTBeCHHO CJIO7KHbIC HOBCPXHOCTH 


BaarogqapHocrH: aBTopbl BbIpaxkaroT OaroqapHoctb B.M. JlaBpiqoBy, I.T.H., upodeccopy, 3aBeyyrolleMy Kadezpor 
«THUG» (®BPBOY BO TOTY, r. XaOaposck) 3a 3HauHMble 3aMe4aHHA WM BaXKHbIC COBCTLI Ip MpoBeyqeHHu 


HCCC HOBaHuA UV odopmMiieHun CTaTbH. 
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Introduction. The reliability of machine parts is determined by such performance properties (PP) of surfaces as 
wear resistance, tightness, strength, quality of coatings [1]. These PP depend on the physico-mechanical and geometric 
parameters of functional surfaces, including roughness [2—4]. 

The analysis of the scientific literature suggests a growing interest in the topic of providing the necessary roughness 
parameters due to the reasonable selection of trajectories of shaping movements and orientation of the sphero- 
cylindrical tool when milling spatially complex surfaces (SCS) [5—7]. Examples of such parts are forming elements of 
die tooling, master models for casting, executive surfaces of gearing [8—10]. 

A number of authors studied the influence of strategies under the milling of SCS and methods of optimizing 
machining [10-12]. However, knowledge about the formation of trajectories does not take into account the 
recommendations for providing the quality of microgeometry of the surfaces of the part. It should also be noted that 
when creating programs for CNC equipment in CAM systems, the limitations of cutting modes are determined 
exclusively using a geometric approach [13, 14]. It does not take into account the influence of the orientation angles of 
the sphero-cylindrical tool relative to the normal plane on the quality of surface treatment, namely on roughness. The 
method of selecting the angles of tool orientation based on empirical models can overcome these disadvantages. Its 
advantages: 

— influence of the tool orientation angles on the surface roughness is taken into account; 

— ability to reasonably select processing paths is supported. 

The study was aimed at the creation of a methodology for selecting the limiting values of the orientation angles of a 
sphero-cylindrical tool to optimize the process of machining spatially complex surfaces. The tasks included achieving 
the minimum values of the amplitude roughness parameter Rz and determining the effectiveness of various machining 
trajectories. 

Materials and Methods. Thus, CAM systems provide forming multi-coordinate machining trajectories with tracking 


of additional parameters, such as collisions, the point of contact between the tool and the part, etc. The sphero- 
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cylindrical tool touches the part at point P; (x;. y;. 2) = Pa (xu. ya. Za). At the same time, it is required to avoid machining 
with the center of the cutter and orient the tool with an angle of inclination of at least 5°-15°. 

In the final operations, the effective cutting speed is determined by the effective diameter. At an equal rotational 
speed, it grows with the increase in the angle of inclination of the tool to the workpiece. An increase in the cutting speed 
generally causes a decrease in the microhardness of the surface, and with an increase in V>75 m/min, the 
microhardness parameters change slightly [12]. The dissipation rate strongly depends on the cutting speed and the 
volume of the material being removed; therefore, cutting-tool lubricant (CTL) is needed to intensify the cutting 
process [15]. 

For the experiments, technological equipment with CNC was used, a five-axis machining center DMU 50 ecoline 
with a maximum spindle frequency of 8,000 rpm. The surface roughness was measured by a Surfcom 1800 D 
profilometer. Sandvik end mills of the R216 series were used for processing 12X18HI10T steel. The material was hard 
alloy 1620 with PVD coating (the closest domestic analogue is T15K6). The diameter was 8 mm, the number of 
teeth — 2. To provide a uniform allowance (a, = 0.2 mm), mechanical treatment with sphero-cylindrical cutters was 
carried out before the final milling operation. 

Research Results. Before determining the angles of inclination, it was required to establish how variable factors 
affected the response function. In this case, we are talking about the surface roughness according to parameter Rz (um). 
To find empirical mathematical models of milling with a sphero-cylindrical tool, we took the independent variables: 
X,— diameter (D, mm) and X2 — the angle of the tool inclination (y, °). The initial data for the analysis were 


considered in previous studies (when applied to tooth £ = 0.4 mm/tooth) [16-19] (Table 1). 


Table 1 
Roughness parameters Rz depending on the angle of inclination and diameter of the tool 
Tool diameter, mm 
Angle, ° 

6 8 10 12 
10 9.33 7.66 5.99 4.33 
20 8.59 7.06 5.53 4.01 
30 7.85 6.46 5.07 3.69 
40 7.11 5.86 4.61 3.37 
50 6.37 5.26 4.15 3.05 


Based on theoretical data on significant factors affecting roughness, linear (1.1) and hyperbolic (1.2) models were 


adopted: 


Rz=Y =atb,X,+b,X,, (1.1) 
hb 

Rz=Y'=a'+b,X,+—. (1.2) 
x 


2 
Here are the calculated natural regression coefficients: a = 13.37; a’ = 10.25; b, =-0.66; b, = —0.58; b, = 0.51. 
The parameters of the two-factor regression equation were estimated using the standard least squares method; 
therefore, for simplicity of presentation, we omitted the formulas indicating the coefficients. Standardized 
B-coefficients: B; =—0.79; B, =—0.58; B, =—0.51. Comparison of the modules of the values of standardized regression 


coefficients B allowed us to conclude that factor X; (tool diameter) had more effect on the formation of the roughness 
parameter Rz than X2 eae angle). Coefficients of multiple, partial, paired correlation and multiple determination: 
Ry x =0.98; tee ~~ 098, he g = 097) tee yp =- 0.95; 


Ry = 0.95; 7, =-— 0.93; r, = 0.85; r , = 0.79; 


° 3 3 
y'xX Y'X)-X> Y'X'5-X, X/X5-Y' 


—_— Fy, =~ 0.585 ry x, = 0.00; ry, = 0.515 


7 YX), 
Ty x =0.00; R2(Y) =0.95; R2(Y') =0.90. 
Comparing the coefficients, we drew the following conclusions. 


When factor X> was fixed at a constant level, factor _X1 most strongly affected (|0.98] > |0.79|). When comparing the 
coefficients of the hyperbolic model, (|0.93| > |0.79]). 
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When factor X; was fixed, the effect of factor X2 on Rz increased for both models: linear |0.97| > |0.58|, hyperbolic 
[0.93] > 0.51]. 

To ensure the uniformity of the microrelief of the surface, the dependence of the feed and the effective diameter of 
the tool (Deap) was established, which varied depending on the angle of processing. We defined the limitations 
associated with the angles of the tool inclination when processing the SCS. To do this, the surface of the part was to be 
divided into sections and the normal angles calculated. If z= f(x. y), then, in general, the orientation of the tool to the 
surface was set by selecting the direction of the normal. 


At cos y=1/|M: 
N= 25 : (2.1) 
ox Oy 
At cos y =—1 /|M: 
N= rer ; (2.2) 
dx Oy 
To determine the inclination angle of the tangent plane, the following equation was used: 
) 2 
tana = |grad(z),| = (= + = ; (3) 
6x by 


where a = |90° — y|. 
AV Nikitenko [20] presented a model for optimizing the orientation angle of a part with corrective angles of 


inclination 4 and B relative to X and Y axes: 


2 2. 
tan a'= [Stan a) + Laer ‘ (4) 
6x éy 


For a special case (Fig. 1), determination of angle A to normal N: 


i = arctan a (5) 
Ax 


| — 
x x X 
i+] 


Fig. 1. Scheme for calculating the angle of the normal: N — normal; 4 — angle to the normal; Ax — selected step along X axis to 
calculate the lengths of the polyline curve segments, mm; Az — distance along Z axis, depending on the step along X axis, mm 


With a discretely defined surface profile, the length of the curve describing the profile geometry: 
S25 AS. (6) 
i=l 


Here, the length of the polyline section AS, = JAx? + Az?. 
The roughness of Rz was considered as an output parameter (Fig. 2), taking into account the limitations associated 
with the trajectories of motion and the angles of inclination of the sphero-cylindrical tool. 
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Fig. 2. Profilograms of surfaces obtained with different shaping trajectories at [ = 35°-45°: 
a — passing milling; b — counter milling; c — top-bottom milling; d — bottom-up milling 
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x, mm 


Top-bottom milling is characterized by the greatest amplitude, the unevenness of the resulting surface profile, and is 


not recommended for shaping with a tool tilted at an angle of 5°-35°. 
The roughness selection map (Fig. 3) for Rz parameter is based on the results of the given and previous 


studies [16-19]. 


Rz, wm 
Gainer 6.1—-10.0 
4.1-6.0 
Top-bottom 
3.0-4.0 
Bottom-up 
No data 


= 10° 20° 30° 40° 50° 60° 70° 80° 
Fig. 3. Roughness selection map 


When using CTL, a film was formed on the contact surfaces of the tool and the workpiece material, which helped to 
reduce adhesive wear. At the cutting speed V> 70 m/min, the effect of dynamic friction was reduced. At the same time, 
the duration of the physico-chemical effect of the medium on the contact surfaces went down, which limited the effect 


of the use of the CTL (Fig. 4). 
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17.5 3155 57.5 77.5 V, m/min 


Fig. 4. Dependence of roughness parameter Rz on the machining speed: a — without CTL; b — with the use of CTL 


The considered technique was aimed not at establishing critical values of possible orientation angles of a sphero- 
cylindrical tool for a specific object, but at achieving roughness parameters taking into account the effective cutting 
speed, feed, and inclination angles for a wide range of parts with concave-convex and linear sections. This approach 
could provide generalizing and clarifying the ways of optimizing machining. In addition to roughness, the limitations of 
the minimum effective cutting speed, depending on the effective diameter of the tool, were analyzed. At the same time, 
the minimum recommended effective cutting speed was (Veq»y) — 75 m/min. 

According to the feed and lateral pitch, the angle of orientation of the tool can correspond to positive and negative 
values. When calculating, it was considered modulo. Based on the calculated data (Table 2), the surface profile (Fig. 1) 
was divided into sections. The normal angles were determined, and the trajectories of the shaping movements were 
assigned to provide the required roughness, taking into account the angles of inclination of the tool. 


Table 2 
Design parameters for determining the optimal angle of inclination of the tool 
n Az, mm Ax, mm 1 a be ____ Veap 
yal yes {=2 at y = 5 
1 0.16 3.59 4.59 6.59 8.59 68.0 
2 0.48 10.83 11.83 13.83 15.83 84.4 
3 0.82 18.25 19.25 21.25 23.25 99.8 
4 1.22 26.01 27.01 29.01 31.01 114.2 
5 1.71 34.34 35.34 37.34 39.34 127.2 
6 2.38 43.64 44.64 46.64 48.64 138.6 
7 3.55 54.82 55.82 57.82 59.82 147.5 
8 6.59 69.21 70.21 72.21 74.21 150.7 
9 4.56 61.26 62.26 64.26 66.26 150.1 
10 2.81 0.25 48.39 49.39 51.39 53.39 143.1 
11 1.98 , 38.41 39.41 41.41 43.41 132.7 
12 1.43 29.71 30.71 32.71 34.71 120.3 
13 1.00 21.73 22.73 24.73 26.73 106.5 
14 0.63 14.17 15.17 17.17 19.17 91.5 
15 0.30 6.87 7.87 9.87 11.87 75.5 
16 0.01 0.32 1.32 3.32 5.32 60.2 
17 0.33 7.52 8.52 10.52 12.52 77.0 
Ni 1.03 22.43 23.43 25.43 27.43 107.8 
Ni+1 1.47 30.46 31.46 33.46 35.46 121.5 
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The measured roughness values, taking into account the recommended angles of inclination of the sphero- 
cylindrical tool and the movement trajectory, are minimal with respect to Rz parameter (from 3 to 6 um). At the same 
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time, these values correlate with data from other studies (Fig. 3). This allows us to conclude that the proposed method 
of preparing control information is adequate. 

Discussion and Conclusion. The presented technique of selecting the limit values of the orientation angles of a 
sphero-cylindrical tool can be used to process the SCS with one tool without replacement, taking into account the 
accepted restrictions. The proposed approach makes it possible to determine the optimal values of the orientation angles 
of a sphero-cylindrical tool, allowing for the cutting speed and achieving the minimum possible amplitude parameter 
of roughness Rz. 

We considered the situation for surface areas with a total angle of 5°-50° at feed /z = 0.4 mm/tooth. In this case, 
machining along trajectories in the passing direction, from bottom to top and in the opposite direction allowed for 
roughness in the range of 3-6 um according to Rz parameter. This was less than the maximum values obtained by 
15-30 %. At angles of 10°—40° and the passing processing direction, the minimum values of Rz — 3—4 um were 
recorded. The trajectory of top-bottom movement was not recommended for use in final operations due to the 
significant height of the Rz profile. At the same time, the values of 4.1-6 um for this trajectory were achieved in a 
narrow range of angles — 40°-S0°. 
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Abstract 
Introduction. \n recent years, one of the main trends in the field of testing road structures has become field study of their 


large-scale models at the accelerated load facility (ALF). It can significantly reduce the cost of selecting the most 
economical and durable pavement designs. However, the results obtained on the ALF are often relative, since they 
practically do not correlate with the results of laboratory and field tests on real objects. This study is aimed at a 
comprehensive investigation of the response of a road structure to a dynamic load, the establishment of patterns of fatigue 
failure of asphalt concrete layers during the accelerated testing and full-scale tests on real objects. 

Materials and Methods. During testing, an accelerated load facility was used, located on the territory of the ShanDong 
Transport University. When conducting field tests, a dynamic loading unit with a falling weight FWD Primax 1500 was 
used, which recorded the deflection bowl on the surface of the structure under study. To record the dynamic response in 
the arrangement of the road structure, a complex of strain gauge sensors was used, which made it possible to register both 
compressive stresses and tensile strains in different layers. The results obtained under natural conditions were compared 
to the results obtained on the mathematical FEM model. 

Results. The research results have shown that the thickness of the lower coating layer is the main factor affecting the 
amount of vertical deformation of the pavement, which must be taken into account at the design stage of the pavement 
structure. Thus, with a thickness of the upper layer of the base of 10 cm, the vertical deformation was 100 um, and with 
a thickness of 20 cm — 55 um, provided that the overall strength of the structure was ensured. The number of load 
application cycles on the ALF had a minimal effect on the selected asphalt concrete samples during split tensile tests. 
Discussion and Conclusion. The adequacy of the results obtained in the course of accelerated testing of road structures 
was shown through a comprehensive comparison of numerical simulation data and full-scale tests, and the adequacy of 
the applied calculation methods was validated. The results of the study can be further applied in the road industry to 
develop and improve the regulatory framework for the design of non-rigid pavement under conditions of increased loads 


and heavy traffic. 


Keywords: asphalt concrete pavement, accelerated loading test, multilayered half-space, stress-strain, temperature 


correction, dynamic response model 


Acknowledgements: the authors appreciate the staff of the Road-testing laboratory, ShanDong Transport University, for 
their assistance in carrying out research at the ALF installation, as well as the respected reviewers for the time and effort 


spent on reviewing this article. 


© Mi Guangcong, Tiraturyan AN, Uglova EV, Vorobev AV, 2023 


Mechanics 


241 


http://vestnik-donstu.ru 


242 


Advanced Engineering Research (Rostov-on-Don). 2023 ;23(3):241-256. eISSN 2687-1653 


Funding information. The research was carried out within the framework of a grant from President of the Russian 


Federation on State support for young Russian scientists — Candidates of science (application MK—242.2022.4). 


For citation. Guangcong Ni, Tiraturyan AN, Uglova EV, Vorobev AV. Study on Dynamic Response Characteristics of 
Different Asphalt Pavement Structures Based on ALF Test. Advanced Engineering Research (Rostov-on-Don). 
2023;23(3):241—256. https://doi.org/10.23947/2687-1653-2023-23-3-241-256 


Hayunaa cmamboa 


MccaeqoBpanne xapakTepHCTHK JMHAMHNYeCKOLO OTKJIHKa JOPOKHBIX KOHCTpyKUMi 


IIpH yCKOpeHHOM TecTHpoBaHHu 


Hal yannynr'®, A.H. Tuparypar’® px, E.B. YraonaO, A.B. BopoGnes*> 
' [lanpayucKuit TpancnoprHsiit yHusepcurer, r. [sunans, nposunuua [anpayn, Kuraiicxaa Hapoguas Pecny6nuKxa 


? Tlonckoit rocyapcTBeHHbIii TexHuueckuit yHuBepcuter, r. Pocros-Ha-Jjony, Poccuiicxas Deyepanusa 


DX] tiraturjan@list.ru 


AHHOTalna 

Beedenue. OHM U3 TiaBHbIX TPCHAOB B OOMACTH MCHbITAHHM FOPOXKHBIX KOHCTPYKUMM B MOCeqHHe TOMI cTasM 
HaTypHble HUCCieqOBaHHA UX KpyMHOMacIITaOHBIxX MOjeselH Ha yCTaHOBKax ycKopeHHoro TecTupoBaHna (ALF). Sto 
TO3BOJIACT 3HAYHTeIbHO YMCHbIUMTb 3aTpaTbI Ha BbIOOp HavOosee IKOHOMMYHBIX HW JOUTOBCYHbIX KOHCTpyKUHi 
HOPOXKHBIX OFeE%*KZ. OFHAKO pe3yIbTaTbl, MOUyYeHHbIe Ha ycTaHOBKax ALF, 3auacTy!0 ABJIAIOTCA OTHOCHTECJIBHBIMH, Tak 
Kak IIpaKTHYeCKH He YBA3bIBAIOTCA C pe3yIbTaTaMM MaOopaTOpHbIX UM MOJICBLIX MCIbITaHHM Ha peaJIbHbIX OOBEKTAX. 
TlosToMy WeuIbiO WaHHOTO MCCIeOBaHHA ABWJIOCb KOMIIJIEKCHOe VW3Y4eHHe OTKIIMKa JOpoOxKHOM KOHCTpyKWMU Ha 
TMHaMV4ecky!0 Harpy3Ky, YCTaHOBJIeHHe 3AKOHOMEPHOCTeli YCTAIOCTHOFO pa3pylWIeHHA ac(asIbTOOeTOHHBIX CJI0eB IPH 
MCIIbITAHHAX Ha YCKOPCHHOe TeCTHpOBaHNe HM Ip HaTyPHBbIX UCHbITAHHAX Ha PeaJIbHbIX OObEKTAX. 

Mamepuanei u memoodvi. Tipu upopeyqenuu UCibITaHHM MCIHONb30BaIaCb yCTaHOBKa YCKOpeHHOTO TeCTHpOBaHHA, 
Haxogamasca B LaHbayHCKOM TpaHCnopTHOM yuHuBepcutete. IloneBbie HCHbITaHvA NpOXoAWIM Cc NpHMeHeHHeM 
yCTaHOBKH JJMHaMHYeCKOrO Harpy2KeHHA C a jarouyHM rpy3om FWD Primax 1500, koTopas ocyulecTBIIAeT perucTpallHio 
yall Wporu6a Ha MOBepxXHOCTH OOcIexyemMol KoHCTpyKunH. Jia perucTpalMu AMHaMM4yecKOorO OTKIMKa B CTpyKType 
WOpOXHOM KOHCTPyKUMH MCHOUb30BaJICA KOMIVIEKC TCH3OMETPHYeCCKHX MaTYMKOB, MO3BOJIAIOWMX OTMe4aTb Kak 
OKMMAIOWMe HallpwKeHHA, TaK M pacTATMBarOllMe epopMalluu B pa3IMYHBIX coax. Pe3yNbTaTbl, MOJYYeCHHbIe B 
HaTYPHBIX YCIOBHAX, ObIIM COMOCTABJICHbI C pe3yIbTaTaMH, WOJyYCHHBIMH Ha MaTemaTHyeckon MK93-mogen. 
Pe3ynomamoi uccnedoeanua. Pe3yimbTatbl MCCeAOBAaHHA NOKa3asIM, 4YTO TOJILMIMHAa BEPXHEFO CJIOA OCHOBAHHA ABJIACTCA 
OCHOBHBIM (aKTOPOM, BIIMAIONWIMM Ha BeIMYMHY BePTHKaIbHO edopMallMH MOpoxXHOTO MOKPbITHA, KOTOPbIT 
HeOOXOAMMO YHUTHIBaTbh Ha CTA MpOeKTHpOBaHHA KOHCTPYKUMU JOpoxKHON oFex Ab. [Ipu TomyuMHe BepxHero ciIOA 
ocHoBaHHa B 10 cM BepTHKalbHad ZedopmMayuna — 100 MKM, a pH TomujMHe B 20 cm — 55 MKM MIpH ycuIOBHU 
oOecreyeHHOcTH OOMei paBHOMpOuHOCcTH KOHCTpyKuUMH. KommuecTBo WHKIOB IpHIOKeHHA Harpy3KH Ha yCTaHOBKe 
yCKOpeHHOro Harpy2KeHHaA HMeeT MHHHMAaJIbHOe BIIMAHHe Ha OTOOpaHHble OOpa3libI accbasbTOOeTOHa Ip MCHbITAaHHAX 
IIpO4HOCTH Ha packOL. 

O6cystcoenue u 3axinrouenue. IlyteéM KOMIVICKCHOTO COMOCTABJICHHA JAaHHbIX YMCJICHHOO MOJeIMpOBaHHA HW HaTYpHbIX 
YCIbITAHHM WOKa3aHa HX TOKCCTBCHHOCTh pe3yJIbTaTaM, TMOIYYeCHHbIM B XOJe€ YCKOPCHHOTO TeCTHPOBaHHA OPO%KHBIX 
KOHCTpyKUMH, OOOCHOBaHa ajleKBATHOCTb IIPHMCHACMBIX PaCUeCTHBIX MeTOAHK. Pe3yIbTaTbI HCCIeqOBaHHA MOTyT OBITS 
TIPHMeHEHEI B JOPOXKHOM OTpaciIM Ia paspaOoTKH MU COBepLIeCHCTBOBAaHHA HOPMaTHBHOM Oa3bI Ip MpoeKTHpoBaHHH 


HO2KCCTKHX JOPO?KHBIX OC7K] B YCJOBHAX MOBbINICHHbIX Harpy30K H HHTCHCHBHOLO ABYW2KCHHA TpaHciopta. 


KoroueBbie cJ1I0Ba: acdbalbToOeTOHHOe TIOKPbITHe, YCTAaHOBKH YCKOpeHHOTO TeCTHpOBaHHaA, MHOrocsIOMHOe 
TIOJTYMpOCTpaHcTBo, Halips2KeHHO-esbopMupoBaHHoe COCTOAHHE, TeMUepaTypHaAr KOppeKTHpoBKa, MOJICJIb 
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BuaarojapHocru: aBTOpbl CTaATbH BbIpaxKalOT HCKPCHHIO!O OaroapHOCTE COTpyHHKaM OpoOo%KHO-HCHbITAaTeIbHO 
jlaOopaTopun MWanpazyucKoro TpaHCIOpTHOrO YHHBepcHuTeta 3a COeCHCTBHE B BBITOJHeCHHH MCCeOBaHui Ha yCTaHOBKe 
yCKOpeHHOro TeCTHpOBaHHA (ALF), a TaK2Ke YBaXKACMBIM PeCUCH3CHTaM 3a BPCMA UH CHIBI, 3ATPAadeHHble Ha paCCMOTpeHHe 
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Introduction. Fatigue damage is one of the main forms of damage to exploited asphalt concrete pavements [1, 2]. As 
a tule, it is associated with alternating loading from the intense multi-cycle impact of the traffic flow. This type of 
destruction is characteristic of highways all over the world and is one of the challenges to which road scientists are trying 
to respond. One of the priorities of modern research in the field of fatigue failure forecasting is the testing of road 
structures at the accelerated load facilities. A certain contribution to the solution of this problem was made by Chinese 
scientists who analyzed the law of variation in the elasticity modulus of the asphalt concrete layer and its deformation 
during accelerated loading [3, 4]. The mechanism of destruction and fatigue characteristics of asphalt concrete mixture 
during operation were investigated in [5]. Its authors obtained a correction factor that took into account the difference in 
the mechanisms of deformation of the road structure under the influence of a test load on accelerated testing installations 
and during real field tests. 

In general, it should be noted that, despite the sufficiently effective simulation of the passage of single loads from the 
traffic flow, reproduced on ALF, difficulties arise with modeling the load distribution over the width of the carriageway, 
under environmental conditions, soil-geological conditions, etc. Thus, the test results at ALF cannot be directly used to 
assess the fatigue characteristics of the actual road surface. They require additional validation under real conditions. 

The authors [6-8] studied the responses of road structures to the dynamic impact of the test load, compared and 
analyzed their changes. The researchers applied methods combining real-time visual observation, instrumental tests for 
deflection, and fatigue failure resistance during loading, and also carried out work to establish relationships between 
stresses, deformations in the pavement structure and the temperature of various pavement structures. 

It should also be noted that similar studies are being conducted in other countries, besides the PR China, including the 
Russian Federation, European countries, and the USA. The results of active research aimed at comparing laboratory 
loading modes to the loading modes corresponding to real road conditions were considered in [9-11]. At the same time, 
there are no similar comparisons of the results of field and laboratory studies to the results obtained on large-scale models 
in a controlled experiment. The results of experimental studies in the field of predicting fatigue failure of asphalt concrete 
layers were given in [12-14]. They comprehensively considered the issues of determining the empirical coefficients 
required to predict the fatigue failure of asphalt concrete layers according to field visual observations, as well as standard 
methods of testing asphalt concrete for fatigue when applying a four-point load or indirect stretching. However, the issues 
of research on large-scale models were also not addressed in them. A number of works deal directly with large-scale 
modeling, but for roadway paving of different designs [15—17]. The comparison of equal-strength structures to different 
layers and their thicknesses is the most interesting from a research point of view, since it will provide establishing patterns 
in the processes of deformation of asphalt concrete associated with its fatigue failure. 

Materials and Methods. Accelerated Load Facility. The Accelerated Load Facility (ALF) is located on the territory 
of ShanDong Transport University. It is a set of equipment for complex testing of road structures (Fig. 1). At the moment, 


it is one of the modeling methods closest to the real conditions for reproducing a moving transport load. ALF is used to 
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study the characteristics of asphalt concrete coatings under plastic deformations, fatigue cracking and downward cracking. 


On the ALF, the axle load and the speed of its movement can be adjusted, which provides the best of all types of test 


equipment compliance with real deformation conditions. 


Fig. 1. Accelerated Load Facility [8] 


Three typical designs of the roadway paving for the device on the test track were selected to study the fatigue resistance 
of asphalt concrete mixtures while in operation and pavement defects in different regions. The test section of the road is 
12 m long and 4 m wide. The road structure is a working layer of the soil of the roadbed, reinforced with cement, with a 
coating of three different types of asphalt concrete: $1, $2, and S3, respectively (Fig. 2). The calculated load reproduced 
by the ALF system is a semi-axial load transmitted through two wheels — BZZ-100. The reproducible load range is 
80-200 kN, the load step is 20 kN, the time interval between load applications is 9 s. 


4m 4m 4m 


5 cm SMA-13 


5 cm AC-20 10 cm AC-20 15 cm AC-20 
10 cm AC-25 


20 cm AC-25 15 cm AC-25 
Soil treated with cement 6 % 


Soil treated with cement 4 % 


Earth road 


Fig. 2. Road structure at the test road section 


A set of equipment was assembled from soil pressure, deformation and temperature sensors installed parallel and 
perpendicular to the direction of loading to provide real-time monitoring of the dynamic response of the structure during 
the construction of the test section of the road (Fig. 3, 4). In the lower part of the monolithic asphalt concrete layer there 
were four sensors of horizontal and longitudinal deformation (Fig. 5). The data collection channel of the corresponding 
sensor and its location are indicated in Table 1. This arrangement of equipment was validated by world experience in 


monitoring the stress-strain state (SSS) of road structures and was described in [18-19]. 
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Fig. 3. Sensor for registering pressure in the soil and on the surface of the base layers [18] 


244 


Guangcong Ni, et al. Study on Dynamic Response Characteristics of Different Asphalt Pavement Structures Based on ALF Test 


CTL 


Fig. 4. Sensor for recording relative tensile strain [18] 


9m 9 
8m 8 
85 
7m 7 
61 
6m [B34] |6 
H- 93 
5m 5 
4m [Pressure [sensor 4 
HH 09 
3m 3 
_W-15(30) [013 
am § c FL24 2 
co 
See 
heh 
0 0 


Fig. 5. The depth of each sensor 


Table 1 
Sensor number and location 

No. Model Location 

Wil B-134 6m 

W2 013 2.5m 

W3 24 2m 

Ww4 09 3.2m 

W5 Y-15 2.5m 

W6 03 5.5m 

Fl Soil pressure sensor 45m 
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To register the dynamic response on the surface of the road structure, in addition to the ALF, a dynamic loading unit 
with a falling load FWD PRIMAX 1500 was used. This unit is a pulse dynamic meter of surface movements of the 
coating, which provides determining the history of movements on the pavement surface under impact by means of 
installed geophone sensors (Fig. 6). Studies have shown that FWD can be used to determine the elastic modulus of 
pavement layers during testing at the ALF [20]. In the course of the research, measurements with the FWD were carried 


out every time after 70,000 cycles on the ALF. The measurement scheme is shown in Figure 7. 


Fig. 6. Falling Weight Deflectometer (FWD) [20] 


Sl 82 83 


Sensor location 


Fig. 7. Schematic diagram of measurements by the FWD 
Split Tensile Tests. To obtain comparable results within the framework of this study, laboratory tests for fatigue 
durability were carried out at the MTS asphalt concrete testing facility [21-24] (Fig. 8). The tests were performed on 
asphalt concrete samples selected from the accelerated testing site in accordance with the PR China standard — JTG E20- 
2011!. When testing asphalt concrete for splitting, the temperature was 15 °C, the loading speed was 50 mm/min. Each 
type of asphalt concrete mix was subjected to four parallel tests. When tested for fatigue durability under controlled 
voltage conditions, the waveform was half-sinusoid with a frequency of 10 Hz. The complete destruction of the sample 


was taken as the criterion of destruction. Each test was carried out for three stress levels and three samples of asphalt 


concrete mix. 
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' Ministry of Transport of the People's Republic of China. JTG E20-2011 “Test Regulations of Asphalt and Asphalt Mixture for Highway Engineering” 
JTG E20-2011 English PDF (JTGE20-2011). (chinesestandard.net) [S]. 
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Fig. 8. Installation of a material testing system [21] 

Research Results. The temperature condition of the road structure operation at the accelerated testing section. 
Temperature sensors were mounted in the test structures of the roadway paving at a depth of 0.2 and 6 cm from the road 
surface, respectively (Fig. 9). 
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Fig. 9. Temperature variation diagram of asphalt concrete layers at different depths 


Automatic data collection was carried out every 15 minutes. Figure 9 shows the curve of temperature variation at 
different depths of the pavement structure. It can be seen that the temperature here changes together with the change in 
ambient temperature. The temperature at a depth of 2 cm from the coating surface is almost during the entire observation 
period above the air temperature and the temperature at a depth of 6 cm. 

Analysis of temperature, stresses and deformations in the lower layer of asphalt concrete. During the tests at the 
ALF, data on tensile deformations at the lower boundary of the lower asphalt concrete layer were collected during each 
test day. To provide the accuracy and reliability of the deformation data, its registration was carried out using a channel 
with a transmission frequency of 2,000 Hz and continuous recording of the response lasting at least three minutes, with 
the possibility of parallel recording of temperature information. 

The maximum number of cycles before the appearance of fatigue cracks during testing at the ALF was 420,000 load 
applications. At the same time, it was found that temperature had the greatest effect on tensile deformations in the range 


from 140,000 to 240,000 applications, as can be seen from the data presented in Figures 10 and 11. 
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Fig. 10. Strain curve depending on the temperature on the lower face of the lower coating layer and the upper coating layer at 
different loading times: a — structure $1; b — structure $2; c — structure S3 
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Fig. 11. Maximum stress value change curve in the lower layer of the coating depending on the temperature at different loading times 


As can be seen in Figures 10 and 11, the deformation of various layers of the pavement changed depending on 
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temperature at different loading times: in S1 structure, the tensile deformation in the lower part of the upper layer of the 


pavement was 19.3—39.7 ps, and the deformation at the lower boundary of the lower layer of asphalt concrete ranged 
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from 58.3 to 75.6 ps. In S2 structure, the range of deformation variation at the lower boundary of the upper coating layer 
was 42.5-64.8 ps, the deformation range of the lower layer of asphalt concrete was 59.2-115.37 ps. Tensile and 
compressive deformations at the lower boundary of the lower coating layer gradually decreased with increasing loading 
time and strongly depended on temperature. Changes on the lower face of the lower layer of asphalt concrete were not 
obvious. In S3 structure, the deformation range in the lower part of the upper coating layer was 76.2—105.2 ws, the 
deformation range on the lower face of the lower layer of asphalt concrete was 92.8—-186.2 us. Tensile deformation of the 
lower coating layer was mainly independent of temperature. 

The range of changes in compressive stress at the upper boundary of the base layer was 30.2—-48.3 kPa; as the number 
of loads increased, the change in compressive stress had a high correlation with the trend of temperature change. 

Analysis of FWD impact test results. The FWD was used, first of all, to determine the actual values of the elastic 
modulus of the layers of the road structure and their variation under testing by the ALF in accordance with the 
backcalculation methodology [25]. 

A total of six tests were carried out by the FWD. At each test point, measurements were carried out three times at the 
same load. Then, the average deflection value at the central loading point for the last two impacts was used as a 
characteristic displacement. The curves of vertical displacement variation at the point of load application depending on 


the number of cycles of application of the load during accelerated tests are shown in Figure 12. 


a 160 

is] 

& 140 

o 

s 120 

S 

5 100 

Se 80 

of 

as. 

3° 40 

a= 

gs 20 

. « 

3 0 50,000 100,000 150,000 200,000 250,000 300,000 350,000 400,000 

a 

a Number of load application cycles 
—e— 0.5m —— LSM rrr 4.5m cece: 5.5m —— 6.5m 
—— 05m —— 1S moe eee 45m sSeiee> 5.5m —— 6.5m 


Fig. 12. Deviation curve of various pavement structures depending on the load time 


The deflection value in the center of the load application by the FWD gradually grew with increasing number of 
loads (Fig. 12). Depending on the temperature and loading time during the first 100,000 load cycles, the change in the 
deflection value was insignificant. The largest change in the magnitude of vertical displacement was observed in the range 
from 100,000 to 200,000 load applications. After 200,000 load applications, the vertical movement stabilized. Since the 
pavement structure was affected by temperature and load time, the structural strength of the pavement decreased 
accordingly, and the amount of vertical displacement increased. 

To verify the results of determining the elastic modulus of asphalt concrete layers using the backcalculation technique 
based on the results of testing by the FWD, additional laboratory examination of asphalt concrete samples was carried 
out’. As a rule, the actual dynamic modulus of elasticity was reduced to a standard value measured at a temperature of 
2 °C according to formula 1: 

Ey) = Ep 0220-1), (1) 


where £29 — inversely calculated module of the asphalt concrete layer at 20 °C, MPa; E7— inversely calculated elastic 
modulus of the asphalt concrete layer at the actual temperature, MPa; 7 — temperature, °C. 


? AASHTO TP 62-2007 Standard Method of Test for Determining the Percentage of Fracture in Coarse Aggregate [S].(D 5821 — 01 (kashanu.ac.ir)) 
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Figure 13 shows the results of determining the elastic modulus of the asphalt concrete base layer depending on the 
number of load cycles reproduced during the accelerated testing. 
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Fig. 13. Variation of dynamic modulus of elasticity of the lower layer of asphalt concrete depending on the number of loading cycles 


Figure 13 shows that as the number of load applications increased, the module of the lower layer of asphalt concrete 
of various pavement structures tended to decrease. At the initial stage of the test, 0-100,000 load applications (the standard 
total calculated axle load was 0—25.6 mln applications), the elastic modulus of the lower layer of asphalt concrete slowly 
decreased with an increase in the number of load applications (the average decrease was 5.23 %). When the number of 
load applications exceeded 100,000, the rate of reduction of the elastic modulus of the lower layer of asphalt concrete 
increased significantly. The modulus of elasticity of the lower layer of asphalt concrete of all structures decreased by an 
average of 54.87 %. 

Split tensile test analysis. As the results of accelerated testing showed, fatigue cracks began to appear by about 
350,000 load application cycle. When the load reached 420,000 applications, fatigue cracks were observed on the surface 
layer of asphalt concrete with an opening width of up to 2 mm, and the test was stopped. 

When sampling cores at the locations of cracks, it was found that cracks, as a rule, appeared on the surface of the 
coating, and did not penetrate into the layers. Most of them were transverse, passing through the rolling bands of the 
wheel load, and only a small part could be attributed to downward fatigue cracks. According to the results of laboratory 
tests of asphalt concrete samples selected at the accelerated testing sections, the dependence of the split strength on the 
number of cycles of the applied load was established. 

Figure 14 shows that the tensile strength during the splitting of asphalt concrete mixtures under operation, subjected 
to different accelerated loading times, changed little. In general, the effect of accelerated loading on the tensile strength 


during the splitting of the asphalt concrete mixture under operation was insignificant. 
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Fig. 14. Strength limit for indirect stretching of asphalt concrete mixtures under operation: 
a — structure S1; b — structure S2; c — structure S3 
Determination of the estimated number of asphalt concrete cycles to failure. At the end of this study, in 
accordance with the regulatory document of the People's Republic of China JTG D50-2017 “Specifications for Design of 
Highway Asphalt Pavement”? and based on the data got by the authors of this study, using the regression analysis method, 
dependence (2) was obtained for predicting the maximum number of cycles before the onset of fatigue failure. To take 
into account the influence of the thickness of the asphalt concrete layer package using the regression analysis method and 


based on the data obtained during accelerated tests, the dependence for the assignment of the correction factor was 


1 4,51 I 1,42 
N= 0265{ 2 (+) : (2) 
g, E 


Np =O) Nay> (3) 


where Ny — number of loading cycles on the actual structure of the roadway paving, cycles; Nai — number of loading 


determined. It is shown in Figure 15. 


cycles on the ALF, cycles; ¢ — tensile strain ratio on the lower face of the asphalt concrete layer, ume; E — elastic 
modulus of asphalt concrete layer, MPa; C; — correction factor between laboratory values and values obtained at a full- 


scale ALF. 
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Fig. 15. Relationship between the ALF correction factor and the thickness of the asphalt concrete layer 


> Ministry of Transport of the People's Republic of China. JTG DS50-2017 “Specifications for Design of Highway Asphalt Pavement. Ministry of 
Transport” [S]. 2017. (STG D50-2017 English PDF (JTGD50-2017) (chinesestandard.net)) 
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The dependence structure corresponds to the general view obtained by the Asphalt institute [14]: 


b c 
aloe 
where a, b, c — empirical coefficients determined under laboratory or field tests; e — ultimate deformation at which the 
tests are carried out; S — material stiffness parameter. 

Analysis of mechanical response characteristics. In JTG D50-2017 “Specifications for Design of Highway 
Asphalt Pavement”, the parameter of fatigue cracking of the bottom layer of asphalt concrete is one of the calculated 
indicators adopted in the design. Due to the greater rigidity of the bottom layer of asphalt concrete, the greatest tensile 
deformations are concentrated at its lower boundary, which is confirmed by the results given in this article. In the 
framework of a numerical experiment, S3 structure is considered as a key control structure for calculating the tensile 
strain in the lower part of the asphalt concrete pavement. The design scheme of the pavement adopted in the simulation 


is shown in Figure 16. 
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Fig. 16. Mechanical response calculation model according to JTG D50-2017 


The calculation of the structure response was carried out on a finite element model of the road structure, while the 
hypothesis of complete adhesion between all structural layers of the pavement was accepted. The calculation was 
performed for a paired wheel of a motor vehicle with an axle load of 200 kN and a pressure of 1.1 Mpa [25-27]. When 
substituting into the mathematical model, the elastic modules of the pavement layers were used as input data, calculated 
according to the test results by the FWD. The results of comparing the road structure response registered during 


accelerated tests and calculated using the finite element method (FEM) are shown in Figure 17. 
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Fig. 17. Comparison of measured and calculated values of the response to deformation of the bottom layer of S3 coating 
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Judging by the data in Figure 17, due to the influence of ambient temperature and the thickness of the bottom layer of 
asphalt concrete, the calculated value of deformation in the lower part is slightly less than the measured value of 
deformation. However, this can be explained by a number of simplifications adopted in the mathematical model. In 
general, the results of full-scale tests and numerical modeling are quite close and are characterized by similar trends, 
which indicates the adequacy of both the full-scale testing and the applied methods of numerical modeling. 

Discussion and Conclusion. One of the basic results presented in this article and differing from the results of previous 
studies [3—-8, 14, 19, 20] is the conclusion about the prevailing influence of the thickness of the top layer of the base on 
the amount of the vertical elastic deformation recorded on the surface of the road structure. It was validated by a set of 
experimental measurements using a dynamic loading unit with FWD Primax 1500. In similar works, it was noted that the 
combined effect of the thickness of the bottom coating layer and the top base layer had the greatest effect on the amount 
of elastic deflection on the surface. 

Also, within the framework of this study, a comprehensive analysis of the test results at the ALF and numerical 
experiments based on known models for predicting fatigue failures was carried out. As a result, dependences (2) and (3) 
were obtained to implement the transition and provide comparability of this outcome. It should be noted that the 
coefficients obtained in the course of this work for dependence (2) differ from similar values known to us obtained at the 


MnRoad test site [14] (Fig. 18). 


Coefficient MnRoad dependence Dependence obtained by the 
authors 
a 0.314 0.265 
b 3.291 4.510 
c 0.854 1.420 


Fig. 18. MnRoad dependence coefficients and those obtained by the authors 


Differences in coefficients are largely due to differences in the rules for the selection of asphalt concrete mixtures, as 
well as regional testing conditions. It should also be noted that these coefficients somewhat tighten the classical 
performance model, reducing the predicted service life of road clothes, in comparison to the MnRoad model. 

Undoubtedly, these results were obtained for a fairly limited sample of experiments and a set of materials according 
to their physical and mechanical properties. They meet the requirements for asphalt concrete mixtures used in the 
construction of highways in China, and are quite close to the mixtures used in the Russian Federation. At the same time, 
it seems promising to use the results obtained. In the future, investigation on calibration of the dependence will be 
continued taking into account data from test sections of operated highways in the Russian Federation and China. 

Laboratory checkout has shown that the split tensile strength practically does not depend on the number of loading 
cycles of the tested asphalt concrete during accelerated testing. 

The results obtained can be applied for further development and improvement of the legal framework of the road 
industry regulating methods of calculation and prediction of fatigue cracking. 
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Abstract 

Introduction. The task of analyzing the stability of plates and shells under creep conditions is critical for structural 
elements made of materials with the property of aging, which are under the action of long-term loads, since the loss of 
stability can occur abruptly and long before the exhaustion of the strength resource of the material. Currently, the issues 
of joint consideration of geometric nonlinearity and creep in the problems of buckling plates remain poorly studied, 
existing software systems do not provide such calculations. The objective of this work is to develop an algorithm for 
calculating the stability of rectangular plates with initial deflection, which are subjected to loads in the middle plane, 
taking into account geometric nonlinearity and creep. 

Materials and Methods. When obtaining the resolving equations, the geometric and static equations of the theory of 
flexible elastic plates were taken as the basis. Physical equations were derived from the assumption that total strains were 
equal to the sum of elastic strains and creep deformations. Finally, the problem was reduced to a system of two differential 
equations, in which the desired functions were the stress and deflection functions. The resulting system of equations was 
solved numerically using the finite-difference method in combination with the method of successive approximations and 
the Euler method. As the boundary conditions for the stress function, the frame analogy was used, as in the case of a plane 
problem of elasticity theory. 

Results. The solution to the problem for a plate compressed in one direction by a uniformly distributed load has been 
presented. The nature of the growth of displacements at different load rates and initial deflection was studied. It has been 
established that when the vertical displacements reach values comparable to the thickness of the plate, their growth rate 
begins to decay even at a load greater than the long-term critical one. 

Discussion and Conclusion. The results of stability analysis using the developed algorithm show that the growth of plate 
deflection under the considered boundary conditions is limited, stability loss is not observed at any load values not 
exceeding the instantaneous critical one. This indicates the possibility of long-term safe operation of such structures with 


a load less than instant critical one. 


Keywords: stability, creep, plate, geometric nonlinearity, physical nonlinearity, initial imperfections, finite-difference 
method 
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AHHOTalna 


Beedenue. 3ajaua ananv3a yCTOMYMBOCTH MIaCTHH HW OOOMOUeK B YCOBHAX MOU3YYeCTH aKTYAaJIbHa [JIA JIEMCHTOB 
KOHCTpyKIUHii H3 MaTepHasioB, OOaqaloulMx CBOMCTBOM CTapeHHA, HaxXOJAMWIMXCA HO WeviCTBHeM JJIMTeJIbHBIX Harpy30K, 
HOCKOJIbKY HOTepx YCTOMYMBOCTH MOXKET MpOHCXOAMTb pe3KO M 3aONrO WO UCueplaHHaA MpOuHOCTHOTO pecypca 
MatepHana. Bompocbl coBMecTHOrO y4eTa TeOMeTPHYeCKOM HeMHeEMHOCTH MU MOU3y4ecTH B 3aayax BbINyYMBAaHHA 
TUIACTHH B HAaCTOAINee BPEMA OCTAIOTCA CIaOO H3YYeHHBIMH, CYWIECTBYIOMMe IPOrpaMMHBIe KOMIJICKCHI He MO3BOJIAIOT 
BBINOJHHTb Tako pacuér. Llenb1o HacTosMel paOoTEI BLICTyMaeT pa3spaOoTKa ajIropuTMa pacueTa Ha yCTOMUYMBOCTb 
TIPAMOYTOJIBHBIX TIACTHHOK C HavaIbHOM MOTHOBIO, HCIIbITHIBAIOMIMX elicTBHe Harpy30K B CpeJ{MHHOM MIOCKOCTH C 
yu4eTOM TeOMeTPHYeCKOM HeIMHeHHOCTH HW MOU3y4ecTH. 

Mamepuanei u memoovi. pu nonyyenuu pa3pelmaroulMx ypaBHeHHii B OCHOBY MOJO%#KeHbI TeoMeTpHyeckHe 
cTaTHueckve ypaBHeHHaA TeOpuu rHOKHX yUpyrHux WiacTHH. Du3st4eckHe ypaBHeHHA BbIBOATCA 13 MpeANOOKeHHA, 4TO 
HOMHbIe WepopMalHu paBHbI CyMMe yupyrux WedopMalui uv epopMalui Nom3yyecTH. OkoHYaTeIbHO 3ayjaya Obuia 
cBeyeHa K CHCTeMe 13 AByx AuddepeHMasbHbIX ypaBHeHHH, B KOTOPbIX B KaYeCTBE HCKOMBIX (PyHKI[MM BbICTYHaIOT 
(YHKUMA HallpwKeHH u npornOa. PemieHue NouyYeHHOM CHCTeMBI ypaBHeHHM BBITIOIHAIOCb YMCIICHHO C MOMOLIbIO 
MeTOJja KOHCYHBIX pa3HOCTell B COUeTAHHH C MeTOZOM MOCIeAOBaTeIbHBIX IpHOMMKeHHM HU MeTOZOM Diimepa. B 
Ka4eCTBe TpaHHYHbIX YCNOBU Id PYHKUMM HallpsxKeHHM MCHONL3yeTcA PpaMHad aHaJIOrMA, KaK B CJIydae TMIOCKOM 
3aja4H TEOpHH yIpyrocru. 

Pezyismamot ucciedosanua. B pamkax TlocTaBieHHO Wem paspadotaH aIropHTM pacueta H pescTaBeHo pelleHue 
3aaqH JIA TWIaCTHHbI, CAKMMACMOM B OJHOM HallpaBJIeHHH paBHOMepHO paciipeyeneHHo Harpy3Kkou. UccneqoBan 
XapaKTep pocta MepeMeleHH pH pa3IM4HOM BeIMYMHe Harpy3KH HW HayabHOH WorMOu. YcTaHoBeHoO, 4TO pu 
JOCTWKEHHM BeEPTHKAIbHbIMH TepeMeLICHHAMH BCIMYHH, COM3MEPHMBIX C TOJIMIMHOM MWIaCTHHKU, CKOPOCTb HX pocta 
HauyHHaeT 3aTyXaTb Jake IPH Harpy3Ke OobUe JIMTebHOM KpUTHYECKON. 

O6cystcdenue u 3akiiouenue. Pe3ynbTaTbI aHaliu3a yYCTOMYMBOCTH C MCMOb30BaHHeM pa3paOoTaHHorO aJiIroputma 
MOKa3bIBaloT, YTO pocT MporHOa WiacCTHHbl pH PaCCMOTPeHHBIX TpaHW4HbIX YCJIOBHAX OTpaHHyeH, WoTeps 
yCTOM4MBOCTH He HaOOaeTCA Mp JNOObIX 3HAYCHHAX Harpy3KH, He MIpeBOCXOMAWHX MFHOBCHHY!IO KPHTHYECKyIO. ITO 
TOBOPHT 0 BOSMOXKHOCTH AJIMTeIbHOM Oe30nacHol SKCIyaTalMu TaKHX KOHCTpyKUMi pu Harpy3ke MeHee MrHOBeHHOM 
KpHTH4eCKOH. 


Kor0ueBble CJIOBa: yCTOH4MBOCTB, TOJI3YU4CCTh, Wi1acTHHa, rCOMeTpH4eCcKaAa HeJIMHeCMHOCTS, dbu3sn4eckaa HeJIMHCHHOCTS, 
HadaJIBHbIle HCCOBCpIICHCTBa, MCTO, KOHCUHBIX pa3HocTeii 


BuaarojapHocrnu: aBTOPbI BbIpaxKaloT OmarolapHOcTE peqakiHv HU peleH3eHTaM 3a BHHMAaTeCJIBHOe OTHOMICHHE K CTaTbe 
HW yKa3aHHble 3aMedaHhA, KOTOPbIe MO3BOJIMJIM HOBbICHTb Ce KadecTBO. 


Aaa wntTupopanna. A3piepn C.b., UenypHenxo A.C. Bsmyanpanve pAMOYrOJIbHbIX TWiaCTHH TIpH HesMHeMHOH 
mom3yuectu. Advanced Engineering Research (Rostov-on-Don). 2023;23(3):257-268. https://doi.org/10.23947/2687- 
1653-2023-23-3-257-268 


Introduction. Much attention is paid to the stability analysis of thin-walled structures in the form of plates and shells, 
since such structures are widely used in construction and other branches of technology [1-3]. One of the challenges in 
the field of calculating plates and shells is the analysis of their stress-strain state under creep conditions, which is 
confirmed by a significant number of works published recently on this problem in domestic and foreign sources. Thus, 
in [4-8], the issues of buckling under creep of composite thin-walled structures were investigated. In [9], the problem of 
stability of functionally gradient plates was considered, taking into account the dependence of material properties on 
temperature. In [10], stochastic analysis methods were applied to the problem of buckling composite plates. In [11-17], 
the issues of stability of viscoelastic plates and shells under the influence of dynamic and tracking loads were discussed, 
and [18] dealt with plates of medium thickness, taking into account the dependence of material properties on time. 
Mathematical difficulties arising in solving these problems led to the fact that numerous researchers limited themselves 
to linear laws of viscoelastic deformation or considered the case of steady-state creep. The finite element method opens 
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up great possibilities in solving the problems of calculating plates and shells taking into account creep. However, modern 
computational complexes, such as ANSYS, Abaqus, LIRA, etc., contain a limited set of rheological models applicable to 
specific materials in a fixed range of stresses and temperatures. There is a need for alternative calculation methods suitable 
for arbitrary laws of viscoelastic deformation, including nonlinear ones. 

This work was aimed at constructing a system of resolving equations for the problem of buckling rectangular plates 
with nonlinear viscoelastic properties under the action of forces in the middle plane, taking into account large 
displacements, as well as an algorithm for its solution. Note that the problem of the stability of structural elements, taking 
into account creep, cannot be solved as a problem of pure stability. Its solution requires the presence of disturbances in 
the form of initial irregularities. Generally, the initial imperfections are given in the form of the initial loss or eccentricities 
of the application of loads. 


Materials and Methods. Let us consider the calculation method on the example of a plate with an initial deflection 
Wo(x,v), compressed by a distributed load p [kN/m] in the x-axis direction and having a hinge support along the 
contour (Fig. 1). 


Fig. 1. Computational scheme 


In the case under consideration, in the presence of creep, if compared to the theory of elastic flexible plates, the 
difference manifests itself only in the form of physical equations. Total deformations can be presented as the sum of 
deformations of the middle plane (passing into the surface) and bending deformations caused by changes in the curvature 


of the middle surface: 
02w 02w 02w 


Vx =Y°-22 ; (1) 


where €, and ¢, — total linear deformations; y,,— total angular deformations; ¢° and ¢° — linear deformations of the 


middle surface; y° — angular deformations of the middle surface. 


For deformations of the middle surface, the equation of continuity of deformations can be written [19]: 


076° eo Ory? _ | (w+ wy) > @ (w+ Wo ) CO? (w+ Wy ) Ow, . O2W, 07 Wy Q) 
Oy? Ox? Oxy Oxdy Ox? oy? Oxdy Ox? dy? | 
For materials with viscoelastic properties, total deformations can be represented as: 
l ee Sige au 
£.= ml ove, ) +E°56, = rie —Vvo, ) +O 5 = ety ag (3) 


where s*,* ,y;, — creep deformations; E — modulus of elasticity; v — Poisson's ratio; o;, oy, Try — values of stress 


J 


components in the corresponding directions. 
Having expressed the stress components in (3) through deformations, we write down the physical relations in the 


inverse form: 


o,= ; 
l-v 


E E 
(e, +Vvé, —(e% +ve',)); Oo, = al + VE, -(e% +ve));t,, = Tew -y%,). (4) 


The relationship of internal force factors and stresses is determined by integral relations: 


hl2 hl2 hl2 hi2 h/2 Al2 
N,= J 6,dz;N,= | o,dz;S= | t,dz;M,= | o,zdz;M,= | 0,2zdz;H = | t,,zdz, (5) 


—h/2 —h/2 —h/2 —h/2 —h/2 -h/2 


where N, and N, — linear longitudinal forces; S — linear shear forces, M, and M, — linear bending moments; 


H — linear torques; h — plate thickness. 
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Next, we substitute (1) in (4), as well as (4) in (5). As a result, we get: 


(6) 


where D = ea cylindrical rigidity of the plate, N* = z ‘a (e% + ve" )de 
12(1-v?) ? e Been l-v? -i2> * pes 
E ap E al2 
N= s\ + ve }dz, S* = ——~ * dz, 
*  J-y2 a . ) Mi+v) ft 
E Al2 E hi2 E h/2 
M*= e. + ve" )zdz, M* = e* + ve* )zdz, H* =—_——_ | _ y-,zdz 
*  J-y? Lf . :) *  J-v? LS :) ai+v) b 


Values N°, N*,S*,M*,M*, H* have the dimension of internal forces and determine the contribution of creep 


deformations to the redistribution of forces. 
Static equations of the flexible plate theory have the form [19]: 
ON, 
ON, ,_8 = 9-5 4 LG 
ox oy ox oy 
2 2 oM , Oo? 0? (wt 0? (wt 
pe ge yy ey CG) Mi) egg (wi) Ma 
Ox? Oxdy = Oy? : Ox? oy? Oxdy 
Here, gq — normal load on the surface of the plate, which is zero in this problem. 
It is possible to satisfy the first two static equations using the Airy stress functions, introduced by the formulas: 


2 2 2. 
ge oe Og r 
ae 


5) 


(7) 
q=0. 


Oxdy 


After substituting the last three equalities from (6) into the last static equation in (7) and taking into account (8), we 
obtain the first resolving equation: 
C20 02 (w+ Wy) C20 02 (w+ Wy) ) C20 0? (w+ Wy) 


DV‘w=q+q* + ————+ ———- 
one Ox? Oy? oy? Ox? Oxdy  OxOy ) 
2 * 2L7* 02M * 
where q* aM: oo 7 + |, 
Ox? Oxdy = Oy? 
To obtain the second resolving equation, it is required to express the deformation of the median surface from (6): 
N,-vN,+N*-vN- 
69 =— Sc Basia aint ee as jane vN* |; 
Eh Eh\ oy? Ox? ° 
N,-vN_+N*-—vN* 2 2 
g0 =—+__+*__*___* = OY a cage ane : (10) 
: Eh Eh\ ox? Oy? : 
2(1+ 2(1l+v 20 
go 2 ogg sO, | 
Eh Eh OxOy 
Substituting (10) into (2), we get: 
2 2 
| up = 6? (w+w) Ow, (wt my) (wt my) Omg OMG (tv) os + 
Eh Oxdy OxOy Ox? oy? ox? Oy? Eh Oxdy (11) 
O2N* O?N)\ @2Nt ea 
Vv x4 ; : 
Ox? oy? oy? Ox? 


Thus, for the problem under consideration, a system of resolving equations is obtained from two fourth-order 
differential equations (9) and (11). Equations (9) and (11) are nonlinear. In the resulting equations, values ® and w are 
functions of the coordinates x, y, and time ¢. Explicitly, there is no time in these equations, the time dependence is laid 
down in creep deformations ¢:, €;,, y%,, which are taken into account by the introduction of integral quantities 


Nt, Nt, S*, Mt, M*, He. 
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For the calculation scheme shown in Figure 1, the boundary conditions are written as: 


2 20 
atx=0,x =a: N, eo pS e 0;w=0;M, =0; 
“Oy? Oxdy ; 
oP Foti) (12) 
aty=0,y=b:N, = 0;S 0;w=0;M,, =0. 
Ox? Oxdy ; 


Equation (11) for small displacements, in the case of a plate made of elastic material, is a biharmonic equation that is 
used to solve the plane problem of the theory of elasticity in stresses. Frame analogy can serve as the boundary conditions 
for the stress function for the biharmonic equation. The plate contour is considered as a frame, and the stress function on 
the contour are equal to the bending moment / in it, and its derivative along the normal to the contour is the longitudinal 
force N. Plots M and N in the frame can be constructed in one of the basic systems of the force method (BSFM). The 
basic system, as well as the diagrams of the bending moment and longitudinal force in the frame are shown in Figure 2. 


Fig. 2. BSFM and diagrams of the bending moment and longitudinal forces 


If the vertical displacements do not exceed a quarter of the plate thickness, then it is possible to assume the forces in 
the middle surface independent of the coordinates x and y (N, =—p,N, =S =0) and use the linearized equation for 


calculations: 


(13) 


The analytical solution to the system of equations (9) and (11) is associated with great difficulties. The authors propose 
to solve this system numerically. The finite-difference method (FDM) was used in combination with the method of 
successive approximations. The Euler method was applied to determine creep deformation in the time domain. As the 
first stage, the solution for the elastic plate was performed. Load p was applied stepwise with small steps. At the initial 
load values, deflections w: were calculated by solving the simplified equation (13). Then, the calculated values w: were 
substituted into the differential equation (11). This provided determining the stress function. The next step was to solve 


the differential equation (9) using the known values of function ®, which made it possible to determine the nodal values 


of deflections w; . After that, the values (11 w, = (w, +wy ) /2 were substituted into equation (11). The iterative process 
was repeated at each step until the relative discrepancy with the norms of the vectors of the nodal values of deflections 
w, and w! was greater than the specified value (the authors assumed it to be equal to 0.1 %). For the second load step, 
the initial value w, in each node was the final result obtained in the first step. The calculation technique in the time 


domain, taking into account creep, was similar. The time interval, at which the process was investigated, was divided into 
steps At. In the case of setting the law of viscoelastic deformation in differential form, the values of creep deformations 


at step ¢ + At were calculated on the basis of the known rate of their growth at time ¢ using the Euler approximation: 


eo. =ert At. (14) 


The block diagram of the creep calculation algorithm is shown in Figure 3. 
Note that the system of equations (9) and (11) provides using schemes of a higher order of accuracy, e.g., the fourth- 


order Runge-Kutta method. At the same time, to achieve the same accuracy of the results, you can set noticeably large 
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time steps. However, with an increase in the step, there is a chance of not catching the effects of unsteady creep at the 
initial moments of time. And with the same time step, the Runge-Kutta method, in comparison to the Euler method, 


requires four times more operations to be performed. 


Solution to the elastic problem 
i=]t=0;e =e) = Vy =0;0=0,,, 


Output of results 


Fig. 3. Block diagram of the creep calculation algorithm 
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Research Results. A polyvinyl chloride polymer plate with dimensions a =2 m, b= 2 m, h=1 cm at E = 1480 MPa, 
v = 0.3 was considered. The nonlinear Maxwell-Gurevich equation was adopted as the law determining the rate of creep 


strain growth: 


0&;, _ iy a _—— 
2E=X,V, J=X); 
oto 
os 6, )- £€;5 
Si =5(5) -60 Wr 2b ip (15) 
* * is * 3 * 
aa] =No ol * ST iis = (o,, 0, ) E,£,, 2 
m 2 max 
where 6, — Kronecker symbol; 5,=J,/3, J,=0,+0, — first invariant of the stress tensor; nj, E. and 


m* — rheological parameters of the material, called the initial relaxation viscosity, the modulus of high elasticity, and the 


modulus of velocity. 
Indices rr in formula (15) indicate the direction of the primary stresses. 
For PVC, the authors took the values of rheological parameters from [20]: E.=5.99-10? MPa, m* = 12.6 MPa, 


ni, =5.44-10’ MPa:s. The shape of the initial deflection wo(x,v) was taken in accordance with the first form of the loss of 


stability of the plate made of elastic material: 
Wy (x9) =fysinsin =. (16) 
a 


For a plate made of elastic material without initial imperfections, the critical load, in the case of the whole ratio of 
sides a/b, was determined from the formula [19]: 
4n?D 
be 
To verify the developed calculation algorithm, the first step was to solve the elastic test problem and compare the 
results to the calculation in the finite element LIRA-CAD package (Fig. 4). The value of the arrow of the initial 
deflection fo was set to 0.15 mm. The grid size when using the FDM was 20x20, the number of load steps was 200. 
When calculating in the LIRA-CAD PC, the plate was divided by triangular finite elements with a triangulation step 
of 0.1 m. The load step size was assumed to be the same as when using the FDM. The value of the critical load for the 
elastic plate, calculated from formula (17), was 1340 N/m. Table 1 shows a comparison of vertical displacements in 
the center of the plate for different load values obtained by the author's method and using the finite element method 
(FEM). The deflections calculated using two alternative methods are quite close, except for the load of 1330 N/m. The 
deviation at this load value can be explained by the fact that when approaching the critical load, the movements rush to 


(17) 


Pr = 


infinity. 


—— oD —>*> EI 6@w6@6S Oa ee 
0 0.0365 0.0456 0.911 1.37 1.82 2.28 2.73 3.19 3.65 


Nonlinear loading 1 
Isofields of displacements by Z(G) 
Units — mm 


= 
YelsX 


Fig. 4. Isofield of vertical displacements in LIRA-CAD PC (p = 1330 N/m) 
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Table 1 
Comparison of calculation results using the author's method and with the help of FEM 
w:103, mm 
p, N/m 
LIRA-CAD the authors’ method 

133 16 16 

266 37 37 

399 63 63 

532 98 99 

665 146 148 

798 218 221 

931 336 342 
1,064 562 578 
1,197 1,163 1229 

1,330 3,646 4,229 


In [21], the possibility of transition from the solution of the elastic problem of calculating plates to the solution at the 
end of the creep process is shown. The value of the long-term critical load p.. can be obtained by replacing the cylindrical 
stiffness D of the elastic plate with the long-term cylindrical stiffness D., which is determined from the formula: 

ahs 


(18) 


1 1 vol 
where a =—+—;B=—+ : 
EE, E 2E, 

For viscoelastic rods and round plates, it has been previously established that in case p<p., the growth of 
displacements in time slows down, and the deflection arrow comes to a finite value. If p =p., deflections grow at a 
constant rate. At p > p., the rate of deflection growth increases. 

The authors also analyzed the nature of the growth of deflections over time for p < po, p = pw and p > p~ for different 
values of the maximum initial loss fo. The deflection curves over time in the center of the plate at p = 0.9 px, p = p» and 


p = 1.1 p»are shown respectively in Figures 5—7. 


w, mm 


1.4 ; T Tr T T T T 7 T 


— f=0.025 mm 


1.2 —— f=0.05 mm 


— f=0.1 mm 


1.0 j—— fo=0.15 mm 
0.8 
0.6 


0.4 


0.2 


Fig. 5. Deflection curves over time in the center of the plate at p = 0.9p« 
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—_ f=0.025 mm 
— f=0.05 mm 
—— fo=0.1 mm 


—— fo=0.15 mm 


0 2 4 6 8 10 12 14 16 18 20 ¢, hour 


— fo=0.15 mm 
—— fo=0.1 mm 
—— f=0.05 mm 
—— fo=0.025 mm 


0 1 2 3 4 5 6 a] 8 9 10 ¢, hour 


Fig. 7. Deflection curves over time in the center of the plate at p = 1.1 po 


Discussion and Conclusion. It can be seen from Fig. 5 that at p < px, the deflection arrow always comes to the final 
value, regardless of the values of the initial imperfections. At the same time, at p > p., the pattern of deflection growth 
obtained in [21] occurs only with small displacements. When deflections reach values exceeding about a quarter of the 
plate thickness, the rate of deformation growth begins to decrease even in the case of loads exceeding the long-term 
critical one. It should also be noted that there is a complete absence of a section with an increasing rate of displacement 
growth for plates with large initial curvatures. The identified effects can be explained by the redistribution of efforts 
N,, N,,S in the middle surface. 

Summarizing the above, we can conclude that the vertical movements of the plate pivotally supported along the 
contour, under the action of a compressive load on one axis, always come to a final value if the load does not exceed the 
instantaneous critical one. In other words, with the considered fastening and loading, the plate is in stable equilibrium 


under creep conditions. 


The obtained equations and the calculation algorithm make it possible to calculate plates made of arbitrary viscoelastic 


materials for any fastening options. The law of the relationship between stresses and creep deformations can also be set 


arbitrarily. 
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Abstract 

Introduction. The performance and reliability of gerotor hydraulic machines depend on the geometry of the cycloidal 
gearing profile. The existing methods of calculating and optimizing the profile parameters are cumbersome, multi- 
criteria and difficult for practical application. Therefore, the problem of creating the methodology for calculating the 
parameters of the gerotor machine profile suitable for engineering calculations at the stage of conceptual design is a 
challenge. In this regard, the objective of this work was to modernize the methodology for designing the geometry of 
the profiles of the end section of hypocycloidal gears used in gerotor hydraulic machines, and to analyze the 
possibilities of their optimization during preliminary design. In the course of the study, the Mathcad computer 
mathematics system was used; numerical experiments were carried out to study the influence of geometric profile 
parameters on the performance and operability of the gerotor hydraulic machine, based on the data obtained and 
analyzed; recommendations for the design of optimal profiles of the end section of gerotor hydraulic machines were 
developed. 

Materials and Methods. Materials included known methods of profile parameters calculation, based on application of 
classical formulas of hypocycloidal equidistant used for outlining profiles of teeth of working elements of gerotor 
machines. The basic research method was modeling the gerotor machine profile using Mathcad computer mathematics 
system. Calculated data were obtained for the selected ranges of varying parameters, processed through the univariate 
regression analysis. 

Results. An algorithm for analyzing the tooth profile smoothness was developed. Two target parameters were defined: 
the cross-sectional area of the end profile, which affects the productivity; the smallest reduced radius of contact that 
determines operability of the working body. A technique for calculating the target parameters at the early stage of 
design was proposed. A number of optimal values of profile parameters according to the criteria of productivity and 
operability of gerotor machine was obtained. The dependences providing the optimum values of profile parameters at 
the stage of designing were constructed. 

Discussion and Conclusion. The developed methodology makes it possible to obtain an assessment of the performance 
and operability of the gerotor hydraulic machine at the design stage of the working body. The research results can be 
used in mechanical engineering when designing gerotor hydraulic machines in order to improve their technical and 
operational characteristics. 


Keywords: hydraulic machine, gerotor gearing, hypocycloid, rotor, stator, contact, open area 
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AHHOTauHA 

Beedenue. (pousBoquTesbHOCTh MH HayjexKHOCTb paOoTb TepoOTOPHbIX THAPOMAIUMH 3aBHCHT OT TeOoMeTpH4ecKHx 
TlapaMeTPoOB IIpoPuIA WHKIOMAaIbHOrO 3allereHuA. CyWjecTByIOWHe MCTOAMKH pacuetTa HW ONTHMH3alMM MapaMeTpoB 
TIpousa TPOMO3AKHe. MHOTOKPHTepHasIbHble. COXKHbIC IA MpakTH4Yeckoro mpuMeHeHua. Tlostomy akTyaJIbHOl 
TIpeACTaBIAeTCA MpoOsIema CO3aHHA METOAMKM pacueta MapaMeTposB Mpodusa padouero opraHa repoTOpHOH MalHHBI, 
IIpHrOAHOH jit MHDKeHEPHBIX pacdeTOB Ha 9Talle 3CKH3HOTO MpoeKTHpoBaHHA. B cBA3H C 9THM IeIbIO JaHHOM padoTsl 
ABIIACTCA MOJCPHH3alHA MCTOAMKUM MpOeKTHpOBaHHA TeOMeTPHH MpoduNel TopleBoro ceYeHuA PHMOUMKMONAasIbHBIX 
3yOUaTBIX 3al[eIVICHHH, UCIONb3yeMbIX B TepOTOPHbIX THUApaBIM4eCKHX MalIMHax. HM aHasIM3 BO3MO%KHOCTeH UX 
ONTHMY3alHH pH IMpesBapHTeIbHOM IpoeKTMpoBaHuu. B xoge UccreqoBaHua Obla HCHOb30BaHa CHCTeMa 
KOMIIbIOTepHOH MaTemMaTHKH Mathcad, mpoBeyeHbI YHCeHHbIe 93KCIepHMeHTbI Id W3yYeHHA BIIMAHHA 
TeOMETPHYCCKHX MapaMeTpoB UpodusieH Ha MIpOH3BOAUTeIBHOCTS UM padoTOcHOcOOHOcTh repoTopHou rHzpaBsM4eckon 
MallIMHbI. Ha OCHOBE MIOJIYYCHHBIX HM MpOaHaIM3HpOBaHHbIX JaHHbIX BbIPAOOTAHbI PeKOMeHAAHH MO IpOeKTHPOBaHHtoO 
ONTHMAJIBHBIX Upodusei TOpueBoro CeYeHHA FTepOTOPHbIX THApPaBJIMYeCKHX MaLIHH. 

Mamepuanot u memoovi. Matepuans BkKm0uaIOT B ceOA M3BECTHbI€ MeTOAHKH pacueTa MapaMeTpoB mpodusa, 
OCHOBaHHBIe Ha IIpHMeHeHHM KaccM4ecKHX (POPMYJI SKBUAMCTAHT THMOWUAKIONA, UCHOb3YeCMBIX JIA OYCPYHBaHHA 
mpoduset 3yObep padOounx opraHoB TepOTOpHbIX MalHH. OcHoBHOM MeTOA UccieqoBaHHa — MOJ{eMpoBaHve 
TIpodusia repOTOPHOH MaLIMHbI C TOMOLIbIO CHCTeMbI KOMIIbIOTepHOM MaTemaTuKH Mathcad. lomyyensi pacueTupre 
JaHHbIe Id BbIOPaHHbIX Wuala30HOB BapbUpyeMbIxX TapamMeTpoB. oOOpaboTaHHbIe MeTOJOM perpeccHoHHoro 
OAHOaKTOpHOroO aHasn3a. 

Pe3yibmamoti uccredoeanua. PaspadotaH asIropuTM aHasiv3a rilaqKOCTH Mpoduse 3yObes. OnpeveseHb! Ba WeseBbIX 
TapaMetpa: IWi0lWayqb CeyeHHa TOpWeBoro Mpo@uNA. BIMAOWat Ha MPOW3BOAMTeIbHOCTb. WM HaMMeHbIUi 
IIpPHBeeCHHbIM payuyc KOHTaKTa. OMpeyeAIOWIMu padoTocnocobHocTs padoyero opraua. IIpeqnoxeHa MeToquKa 110 
pacueTy lJeJIeBbIX HapaMeTpoB Ha paHHel cTaquM upoekxTupoBanua. TlomyyeH pay ONTHMaJIbHbIX 3HayYeHHi 
TapaMeTposB npodusa 10 KPHTepHAM MPOM3BOAHTebHOCTH UM padoTocHocobHocTH repoTopHon MawuuHbl. [locrpoenst 
3aBHCHMOCTH. MOSBOJIAIOWIMe YCTaHOBHTb ONTHMaJIbHble 3Ha4eHHA apaMeTpoB UpopusA Ha CTaqMu MpoeKTHpoBaHHA. 
O6cystcdenue u 3akniouenue. PazpaboTaHHad aBTOpaMH MeTOHKa JlaeT BO3MOXKHOCTS Ha JTalle MpoeKTHpOBaHHA 
paOoyero opraHa MOJyYHTb OWeCHKy MpOu3BOANTebHOCTH HM padoTocMocoOHOcTH TepoTopHol TrHApoMalHHEl. 
Pe3yibTaTbI MCCIeOBAHHM MOTYT ObITb MCMONb30BaHbI B MalIMHOCTpOeHHH TIpH IMpoeKTHpOBaHHH TepOTOPHBIX 
THAPOMallMH JIA yIYaWeHHA HX TeXHMKO-93KCIIyaTalMOHHBIX XapakTepHCTHK. 


KoroueBple CJI0Ba: THApaBlIndeckadA MalliwHa, TepoOTOpHoe 3allemieHve, THNOUNKIONAa, poTOp, CTaTOp, KOHTAKT, 
TJLOMWayb 2KUBOTO CC4YCHHA 


BaarojapHoctrn: aBTOpbI BbIpaxKaroT OlaroapHOcTE BC€M YyU4aCTHHKaM JOrOBOpa 3a TIOMOMIb B TIpoOBeyeHHu 
CTCHJOBBIX HMCIIbITAHHH. 
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Introduction. Gerotor hydraulic machines are increasingly being used in the construction and operation of oil and 
gas wells. The reason is their structural and operational advantages, the primary of which are low delivery unevenness 
for pumps or torque for motors, compact overall dimensions. 

In the scientific literature, a large place is occupied by studies devoted to the operation of gerotor hydraulic 
machines. The strong geometric effects of the rotor profile on the performance and reliability of screw gerotor hydraulic 


Kireyev SO, et al. Optimization of Geometric Characteristics of Cycloidal Profiles of Gerotor Hydraulic Machines 


machines were noted by D. and F. Baldenko [1]. The basic unit of such a machine is a multithread screw gerotor 
mechanism — a cylindrical planetary gear train of internal gearing, consisting of a stator and a rotor, with a difference 
in the number of teeth equal to one. The stator profile in the end section contacts the rotor profile, forming closed 
platforms. The authors also analyzed the technology of manufacturing gears. 

R. Iyer, V. Jatti and others investigated the feasibility of mathematical models in the design of screw pumps [2]. 
They also presented the results of modeling various sizes of twin-screw pumps when pumping water, a mixture of water 
and air. Hong-Seok Park and Xuan-Phuong Dang proposed a CAD model optimization process using the results of CAE 
modeling and Mathcad calculations [3]. Other authors (A. Lebedev. S. Kireev. M. Korchagina. S. Vlaskin. A. Efimov) 
noted in their publications the effectiveness of designing various objects and processes based on parametrization and 
subsequent optimization of the design according to one or more target parameters [4—6]. 

In their work, F.D. Baldenko and Yang Yao have emphasized that the open area does not depend on the relative 
position of the rotor and stator profiles and determines the pump performance. They analyzed the effect of 
dimensionless geometric profile parameters on the open area [7]. 

Jie Chen, He Liu, Fengshan Wang, Guocheng Shi, Gang Cao and Hengan Wu developed a finite element model of 
the working body of a screw pump with a surface drive for the oil and gas industry [8]. Based on the results of modeling 
and experimental tests, they concluded that the gaps in the coupling and the thickness of the stator were the two key 
factors affecting the optimized target parameter — the volumetric efficiency of the pump. 

New data were obtained by J. Gamboa, A. Olivet and S.H. Espin when modeling the working body of a multilobe 
screw pump [9]. These authors determined a significant relationship between the sliding clearance area and the pressure 
drop associated with the mechanical properties of the stator material. 

Yang Yao and F. Baldenko investigated the effect of dimensionless coefficients of the end profile of a gerotor 
hydraulic machine on the characteristics of cycloid engagement, which provided optimization based on geometric and 
kinematic criteria [10]. I-A. Lyagov, F.D. Baldenk, A.V. Lyagov, V.U. Yamaliev and A.A. Lyagova proposed a 
methodology for selecting the optimal configuration of power sections when designing sectional downhole motors, 
noting the urgency of the problem of reducing the overall parameters of the screw pump length [11]. 

Yu.A. Korotaev, A.N. Alpatov, A.V. Sobolev and others have come to the conclusion that the most effective is mesh 
profiling from the initial contour of the rail outlined by the equidistant shortened cycloid [12]. With regard to the 
profiles obtained from the equidistant of an ordinary hypocycloid, the authors pointed out their major drawback: they 
were built as special and did not provide for tension in the engagement. Yu.A. Korotaev and D.A. Goldobin drew 
attention to the problem associated with the need to increase the performance of screw pumps, and the resulting 
unacceptable rigidity conditions under the production of stators and rotors [13]. 

A.F. Minikaev, V.A. Pronin, D.V. Zhignovskaya, Yu.L. Kuznetsov, F.D. Baldenko, A.E. Kovalenok developed a 
methodology for studying the stress-strain state of the stator during contact interaction with the rotor of a screw working 
body of a hydraulic machine using CAE modeling. They have stated that the classical Hertz theory of contact 
interaction cannot be applied to determine contact stresses in the rotor and stator due to a number of features: change in 
the curvature of the contacting surfaces of the rotor and stator; different thickness of the elastic lining along the 
circumference of the housing; presence of preload in a pair; complex planetary motion of the rotor, accompanied by a 
combination of friction rolling and sliding; displacement of the radial force vector relative to the normal at the point of 
contact. 

Summing up the review of scientific publications on the design of a gerotor hydraulic machine, it should be noted 
that most researchers in this field strive to improve the efficiency of this machine by increasing the open area and 
reducing the overall dimensions of the working body. Optimal results were achieved by applying and varying the 
coefficients of off-centroidness, tooth shape, relative sliding speed, rail displacement, etc. Calculation methods were 
cumbersome, they were difficult to apply in practice. End profiles were built, as a rule, on the basis of shortened 
hypocycloids, which caused problems with providing kinematic accuracy of movement, contact interaction of the rotor 
and stator. Some authors (Yu.A. Korotaev, A.N. Alpatov, A.V. Sobolev, et al.) pointed out the impossibility of creating 
tension in the case of using an ordinary hypocycloid, but tension could be created through changing the equidistant 
distance of the rotor and stator profiles. 

Materials and Methods. To obtain profiles of the rotor and stator in the end section of the working body of the 
gerotor hydraulic machine, we specify the following set of initial parameters (Fig. 1): 

1) R,— profile contour radius, which determines the diameter of circle C), describing the profile of the stator (mm). 
In fact, this parameter assigns the radial dimension of the gerotor machine, which is particularly important in the oil and 
gas industry when the projected object needs to be placed, e.g., in the trunk of the production well column; 

2) z, — number of rotor teeth; 
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3) R, — radius of the roller, which determines the distance of the equidistant line of the rotor and stator profiles 


relative to the corresponding hypocycloids (mm); 

4) A — tension in the coupling of the rotor and stator (mm). 

Recall the principle of formation of hypocycloidal engagement profiles, which consists in the fact that inside circle 
C, (base circle of the stator hypocycloid), circle C, (base circle of the rotor hypocycloid) rolls without sliding. The 


center of circle C, is always at distance e (eccentricity) from the center of circle C,,. 


Stator hypocyploid 
“| Mi 


val 


120; 


Stator profile 


-150 -120  -90 -60 -30 0 30 60 90 120. 150 
Fig. 1. Diagram of formation of the hypocycloidal profile of the end section of the working body of the gerotor machine, 
z,=4, R, =150 mm, RF, = 30 mm 
Inside circle C,, circle C, (generating circle of the rotor hypocycloid), one point of which forms the rotor 
hypocycloid, rolls without sliding. Inside circle C,, circle C, (generating circle of the stator hypocycloid), one 
point of which forms the stator hypocycloid, rolls without sliding. Note that the radii of circles C,; and C, are equal 


to eccentricity e. The rotor and stator profiles are formed as equidistant lines deposited at distance R, from the 


corresponding hypocycloids. 
In this setting, the rotor and stator profiles provide an ideal contact under eccentric rotation of the rotor. When 
required to provide tension/clearance in the engagement, the equidistant distance R. for the rotor profile should be 


changed by A. 
Let us move on to the mathematical description of the method of forming the profiles of the end section of the rotor 


and stator. The number of stator teeth is equal to 
2,22, +1, (1) 
The radius of the base circle of the stator hypocycloid R,, is determined from the condition (Fig. 1): 
R,, =R, -R.. (2) 
We determine eccentricity e of the transmission from the condition that the length of circle C, is equal to integer 
z, of the lengths of generating circle C, with radius e: 


2pR,, =2 pez, 


bs 


then, 
e=R,./z,. (3) 
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Note that for z, = 1 (Moineau mechanism) according to other formulas: 
e=R, /3uR, =2e, if z, =1. (4) 


The radius of the base circle of the rotor hypocycloid R, is equal to (Fig. 1): 


br 

R,, =R,, —@. (5) 

We introduce a function to determine the radius of the base circle of the rotor or stator, fulfilling conditions (2) 
and (5): 

R(z)=R, —R, —e-(z, —2Z). (6) 


Parametric formulas of the hypocycloid of the rotor and stator in the Cartesian coordinate system XOY: 


x(@, Z) = R(z)-| (1-1/z)cos(@/z)+(a/z)cos((1-1/z)9) | 


: 7 
y(Q, 2) = R(z)-| (1-1/z)sin(@/z)-(a/z)sin((1-1/z)@) | ” 


where z — number of branches of the hypocycloid of rotor z, or stator z,, a — velocity factor of the hypocycloid. In 
the case of an ordinary hypocycloid in theory, a =1. However, in practical calculations using computer mathematics 
systems (Mathcad. MATHLAB) at a =1, difficulties arise with calculating equidistant profiles of the rotor and stator. 
Therefore, in the model, we assume a = 0.99, which is quite justified to achieve our goals; @ — rotation angle of the 
radius vector drawn from the origin to each point of the hypocycloid. If the hypocycloid must be constructed from 
N,, points, then 

o(i,z)=2-n-i-z/N,. i=0...N,. (8) 


P 


The profiles of the rotor and stator are determined from the well-known equidistant trajectory formulas relative to 


the hypocycloids of the rotor or stator at distance R, taking into account the tension in the coupling A. 


X,,(9,2,A) = X(9,2) +(R, +A): oe (22 2) {ae 2) 


0) 0) (9) 
¥,(9,2,A) = (,2)-(R, +A): “e2 {2 1{Be2) 


During transmission operation, the center of the rotor moves along a circle with radius e and rotates by angle 
a=0...27. At the same time, all points of the rotor make an additional rotation by angle wy, determined by the function 

yw(a) =a-(l—z, /z,). (10) 

The rotor rotation model is described by parametric coordinate transformation functions when moving and rotating 


an object: 


X,(a,9, A) =e-cos(a) + X (9, 2,, A) cos(y(a)) — ¥, (@, z,,A) sin(y(a)) 


Y,(a,0,A) =e-sin(a)-if(z, =1,0,1) +X, (@,2,,A)sin(w(a)) + Y, (@,2,,A) cos(w(a)) ay) 


In expression (11), a conditional function is applied that zeroes the movement along Y axis in the case of rectilinear 
movement of the rotor in the stator cavity at z, =1. 

The experience of using this technique has shown that the equidistant profiles of the rotor and stator with an increase 
in the equidistant distance or radius of the roller R, can lose smoothness. This occurs at the points of transition from the 
equidistant circle to the equidistant hypocycloid. Figure 2 shows the problem of loss of smoothness of the stator profile 
line in the vicinity of +1 mm from the transition point of equidistant lines. 

In mathematics, the formula for the function smoothness criterion given in parametric form is known. For the stator 


profile, such a function has the form 
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OY, (P,2,0) 6X (P20) 


12 
ap ap (12) 


F.(p,2) = 


r) 


where z= z,; 
p — angle of the radius vector of the stator profile point. 
It should be noted that function F. behaves unstable in the entire range @(i,z,), and it is quite difficult to find any 


features of its behavior at the transition point of equidistant lines of interest to us. Therefore, the behavior of function 
F. was studied only in the vicinity of the transition point when the angular parameter changes 


p=t/z, —m...m/ z, +m, (13) 


where m — angle of deviation from the angle of the radius vector of the transition point. 
m=2-p-u, (14) 


where u — part of the angular sector of the stator profile to study smoothness in the vicinity of the transition point. 
Value u (for z, =2..12) is determined from the empirical formula 


u=0.12 —0.01-(z, —1). (15) 


Figure 2 shows that with an increase in the radius of the roller in the vicinity of the transition point of equidistant 
lines, the stator profile turns from smooth to self-intersecting. The behavior of function F, also changes, and when 


the self-intersection effect of the profile appears, a loop occurs (Fig. 2 c). 

This effect is taken as a basis for the stator profile smoothness control algorithm proposed by the authors of the 
article (Fig. 3). It should be noted that, judging by numerous measurements, if the stator smoothness criterion is met, 
then the smoothness of the rotor profile does not need to be checked, since it is built on a hypocycloid with a smaller 


value of the radius of the base circle. 


Yp, Yp, Yp, 
mm mm mm 
Transition point Transition point Transition point 
39.82 49.48 59.13 
39.15 48.81 58.46 
111.1 111.7 112.3, Xp,mm 103.8 104.4. 105.0 Xp, mm 96.67 97.34 98.01 Xp, mm 
F, F, F, 
Smoothness is available Smoothness is still there Smoothness broken 
0.34: 0.34 
3—$$— 21.33, 
4 | 
111.1 111.7 112.3 Xp,mm 103.8 104.4 105.0 Xp,mm 96.67 97.34 98.01 Xp, mm 
a) b) c) 


Fig. 2. Stator profile in the vicinity of +1 mm from the transition point of equidistant lines (upper graph) and values of 
function Fs (lower graph) at z, =3 and R, =150 mm, for R, :a—40mm; b)— 50 mm; c — 60 mm 
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Xp(p,2Z,A),Z,,m,N,, A 


yes no 


Cycle p=1/z, —m,1/z, +m, step Ap 


no 


Smoothness not Smoothness broken 


broken 


Xp(P, 2,0) 2 
Xp(p + Ag, z,,0) 


v yes 
aay 


Fig. 3. Algorithm scheme for monitoring the stator profile smoothness 


Control of the stator profile smoothness should always be carried out when designing the geometry of the end 
section of the gerotor machine. In order to study the parameter of the maximum radius of roller R below which the 


smoothness of the curve was not violated, the authors conducted a series of numerical experiments where various values 
of the contour radius R, and the numbers of rotor teeth z, were specified. 

The first target parameter of this study is the open area of the working body S. This parameter is included in the 
formula for feeding the gerotor machine, and it determines its performance. With the available parametric functions 
describing the profiles of the rotor and stator, the task of determining S is reduced to calculating the difference in the 
areas of the parametric curves of stator S$, and rotor S (9). The formula for the area of the parametric curve is known 


from mathematics; therefore, in the authors' methodology, it looks like this: 


22°20 oY, (pz 0) 2-20 oY, (p,Z ,A) 
S=8, -S. = [0.2.01 - ap ; dp a A, (P.2,,A) P ap Lp (16) 


The second target parameter of the study is the reduced radius p,, of curvature at the contact interaction of the rotor 


and the stator. Recall that according to Hertz theory, when two cylindrical surfaces come into contact, the reduced 
radius is calculated, which is in the denominator in the contact stress formula. Consequently, the magnitude of the 
contact stresses, along with the physical and mechanical parameters of the materials of the contacting surfaces, has the 
largest value at the lowest value of the reduced radius p,,,, determined from the formula: 


Pap =(P, "P.)/(P, +P,)- (17) 


where p,— radius of the equidistant circle of the rotor, equal to parameter R; p,— radius of the stator profile 
curve. 

Note that if p, >1, then the numerator in formula (17) always grows faster than the denominator; therefore, the 
smallest value 9,,, ni, 18 completely determined by the parameter of the smallest radius of the stator p., nin: 

In the end section of the working body of the gerotor machine, at any rotation angle of the rotor, continuous contact 
of its equidistant circle and a complex curve of the stator profile is made. Based on the results of observations of the 
contact interaction dynamics, it can be concluded that the greatest stresses occur at the transition points of the 
equidistant hypocycloids in the equidistant circumference of the stator profile. Consequently, the smallest radius of the 
stator profile line p,, ,,,, 18 also located in the transition zone of equidistant lines and determined by the formula known 


from mathematics: 


3 
OX ,(P.22,0)) | ( 8¥,(p.22.0)) 
op op 
Ponmin (p) = = (18) 
OX , (P22,9) . O7Y, (Pp, 2750) aye O?X ,(P, 21,9) k OY, (Pp, 27,9) 
Op Op? Op? Op 
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where p=x/z, angular parameter of the stator profile point at the junction of equidistant lines. In case of loss of 
smoothness of the stator line, parameter p,,,,,, tends to zero. 

The technique described above was implemented in the Mathcad system of mathematical calculations. When 
performing calculations, the hypocycloid velocity factor a=0.99, N, = 40,000. 


At the first stage of research. it was planned to study the maximum radius of roller R at which the stator profile 


line smoothness criterion G was zero. Variable parameters — contour radius R, (100, 150, 200 mm) and the number of 
rotor teeth z, (from 1 to 7). Calculations were performed with a sequential increase in parameter R. from 5 mm until 


the smoothness criterion G was trigged. The investigated parameter R was fixed with an accuracy of 1 mm. 


r max 


At the second stage, a series of calculations were carried out to study the effect of the roller radius R, (from 5 mm 
to R.,. ), the number of rotor teeth z, (from 2 to 7) and the contour radius R, (100, 150, 200 mm) on the open area S 


and the minimum reduced contact radius ,,, nin- 

The described method of calculating the profile geometry of the end section of the working body of the gerotor 
machine can be implemented in any other programming system or in a computer mathematics system. 

Research Results. The data of numerical experiments to determine the parameter of the maximum radius of roller 


R. ax» at which the stator profile line smoothness is preserved, are shown in Table 1. 
Table 1 
Maximum radius of roller R, ,,,, (mm) at different contour radii R, (mm) 
a R, = 100 R, =150 R, =200 
1 66.6 100 133 
2 45 68 91 
3 27 41 55 
4 20 31 40 
5 15 23 31 
6 12 19 25 
7 10 16 21 


Based on the data in Table 1, graphs were constructed (Fig. 4 a). It can be seen that the estimated parameter R, ,.. 
grows proportionally with increasing contour radius R,. We introduce into consideration the relative parameter 
R, = Rmx /R,, Whose graph is presented in Figure 4 b. 


vr max 


r max ? mm Ro 
140 0.9 
120 0.8 
0.7 1 
100 L 
0.6 
80 ; 0.5 
60 P 0.4 
2 
40 0.3 / D 
20 a4 
0 | 0 
1 3 5 7 71, pes 1 6 11 16 71, pes 
a) b) 


Fig. 4. Dependence of maximum radius of the roller on the number of teeth of the rotor: 
a— R,,,, and R,: 1 — 100 mm; 2 — 150 mm; 3 — 200 mm; b — R_ (1) and trend line (2) with extrapolation to z, =16 


r max 
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The empirical formula from which it is possible to determine the maximum allowable radius of roller R has the form: 


r max ? 


RK, =0.768z,°"" (19) 


At z,=15, R,=150mm, Rk, = 0.0517 mm, R.,,,, = 7.76 mm, the profile smoothness is achieved at R-= 7.5 mm. 


vy max 


At z,=10, R,=200mm, R = 0.0775 mm, R,,,,, = 15.5 mm, the profile smoothness is achieved at R-= 15.5 mm. 


Thus, we can conclude about the satisfactory accuracy of empirical formula (19) for predicting the maximum allowable 
parameter R,- max. 
It should be noted that formula (19) was obtained at the values of the velocity factor a = 0.99. A decrease in a causes 


the smoothness of the profile with a more significant increase in parameter R,. However, the shape of the cycloidal 


engagement tooth changes, and the conditions of contact interaction between the rotor and the stator are violated, 
specifically, in those places where the rotor should touch the stator with the sides of the teeth. Calculations show that at 


a=0.99 and R, = 38 mm, the smoothness criterion is at the limit of feasibility. Though at a=0.90 and R, =38 mm, 


good smoothness of the stator profile is achieved, but the kinematic contact of the parts of the working body of the 
machine is disrupted, which can be restored only by a special selection of the values of the hypocycloid parameters of 


the rotor and stator. Therefore, the study of the influence of the velocity factor on parameter R is beyond the scope 


vy max 


of this research. 


Table 2 presents the data of calculations of the target parameters of the studies for the number of rotor teeth z, =3 


at different values of contour radii. The obtained values of the open area S and the minimum reduced contact radius 


Pupmin 10 the last two columns of Table 2 are given in a dimensionless form for ease of analysis. 
Table 2 
Calculation data of target parameters at z\=3 
R,, mm Pup mins MM S, mm? S/max(S) Pap min » (M0AX( Py min 
R, = 150 mm 
5 4.72 18,038 0.9712 0.2656 
10 8.83 18,256 0.983 0.4969 
15 12.28 18,415 0.9915 0.6911 
20 14.98 18,513 0.9968 0.843 
25 16.85 18,572 1 0.9482 
30 17.77 18,530 0.9977 1 
35 17.63 18,450 0.9934 0.9921 
40 16.28 18,307 0.9857 0.9162 
R, = 200 mm 
10 9.14 32,275 0.9786 0.3857 
20 16.37 32,738 0.9926 0.6909 
30 21.36 32,982 1 0.9015 
35 22.9 32,980 0.9999 0.9664 
40 23.69 32,943 0.9988 0.9998 
45 23.695 32,844 0.9958 1 
50 22.82 32,685 0.991 0.9631 
55 20.98 32,466 0.9844 0.8854 
R, = 100 mm 
5 4.57 8,068 0.9773 0.3857 
10 8.18 8,184 0.9914 0.6903 
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R.,mm Pap mins "1M S, mm? S/max(S) Pap min> /MAX( Pry min ) 
R, = 150 mm 
13 9.83 8,225 0.9964 0.8295 
15 10.68 8,245 0.9988 0.9013 
17 11.32 8.255 1 0.9553 
20 11.85 8,235 0.9976 1 
24 11.64 8,189 0.992 0.9823 


Figure 5 shows graphs of variation in the target parameters of studies at z,=3. It is easily seen on them that with an 
increase in the radius of roller R., there is a slight change in the open area S and a sufficiently strong change in the 
reduced contact radius, whose maximum value falls on a fairly narrow range R.,. This is of the greatest interest to the 
designer, since in this case, the contact stresses are minimal. Calculations for other values z, (from 2 to 7) showed 


approximately the same qualitative picture as in Figure 5. 


1.0 a 7a i ts 1.0 one eFeege 1.0 ee oe 
09 > y 0.9 « a) ie r @ 
0.8 2 ane: le 
fl 0.8 : 0.8 : 
0.7 6 Fa 2 - 
ae ? 0.7 eo 0.7 PZ. 
i ‘ 0.6 0.6 
0.4 : 0.5 0.5 
i 0.4 ; 0.4 : 
03 | ¢ é é 
0.2 0.3 0.30 
0 10 20 30 &,mm 0 20 40 R,mm 0) 10 20 R,mm 
a) b) o) 


Fig. 5. Dependences of the target parameters of the study on the radius of roller R. at z1=3.a=0.99. A=0: 1 — S/max(S); 


2 — Pap min MAX( Pry min ) (Table 2); R, . mm: a — 150; b — 200; ce — 100 mm 


Table 3 shows the calculation data of the optimal radii of rollers R, op for various parameters z, and R,. As a result 
of the analysis, it turned out that parameters R,5,,and R, for one z, change almost the same. Therefore, if you enter the 


parameter of the relative optimal radius of the roller 


R,, =R,,, /R,, (20) 


ro ++ opt 
then, with the help of a single-factor dependence of R,, = f(z,) )-type, it is possible to determine the optimal radius of 
roller R, at the preliminary stage of designing the working body of the gerotor machine. 
As a result of similar reasoning, for parameters S and p,,,,,,,, given in Table 3, dimensionless parameters of the 
relative reduced radius of contact p,,,, and relative open area ,,,,,,, were also introduced (Table 3): 
S,=SI RESP apo =Tapmin / Ry (21) 


New relative values are also presented in Table 3. 


Kireyev SO, et al. Optimization of Geometric Characteristics of Cycloidal Profiles of Gerotor Hydraulic Machines 


Table 3 
Optimal radii &,.,,, and relative engagement parameters 
R,, mm Zi R. op? TAM S, mm? Pup mine 1M R,, Papo S, 
200 2 68 36,103 43.2 0.340 0.216 0.903 
150 2 51 20,318 32.2 0.340 0.215 0.903 
100 2 34 9,015 21.65 0.340 0.217 0.902 
100 3 21 8,235 11.85 0.210 0.119 0.824 
150 3 32 18,530 17.8 0.213 0.119 0.824 
200 3 42 32,800 23.1 0.210 0.116 0.820 
200 4 28 30,015 15.32 0.140 0.077 0.750 
150 4 21 16,880 11.5 0.140 0.077 0.750 
100 4 14 7,495 7.63 0.140 0.076 0.750 
100 5 10 6,860 5.35 0.100 0.054 0.686 
150 5 15 15,441 8.1 0.100 0.054 0.686 
200 5 20 27,452 10.81 0.100 0.054 0.686 
150 7 9 13,050 4.7 0.060 0.031 0.580 
200 7 12 22,100 6.2 0.060 0.031 0.582 
100 7 6 5,780 3.2 0.060 0.032 0.578 


It should be noted that all newly introduced relative values with a high degree of confidence obey approximating 


formulas that include one factor — the number of rotor teeth z,. 


The formula for the relative optimal radius of roller R,, has the form 


R,, = 0.9331-z, 3” (R? = 0.9965), 


where R? — value of the approximation accuracy. 


The formula for the relative reduced contact radius p,,,, has the form 


Papo = 0.635-z, '*(R? = 0.9995). 


The formula of the relative open area S, has the form 


S, =1.0175—0.0641z, (R? = 0.9911). 


(22) 


(23) 


(24) 


Let us consider an example of the application of the results obtained in the design of the geometry of the end section 


profile of the gerotor machine. We specify the following set of initial data: contour radius R, =120 mm, number of 


rotor teeth z, =4, tension in the coupling A=0 mm (empirical formulas are obtained at zero tension). The relative 


optimal radius of roller R,, according to formula (22): 


R,, = 0.9331-4-3” = 0.135. 


Absolute value of the roller radius R, = R,,-R, = 0.135-120 =16.2mm. 


Relative reduced contact radius p,,,, from formula (23): 


Ppp = 0.635-4-"™ = 0.076. 


Absolute value of the minimum reduced contact radius 


p np min 


Relative open area S$, from formula (24): 
S, =1.0175—0.0641-4 = 0.761. 


Absolute value of the open area 


= Tipo 


R, = 9.076 -120 = 9.12 mm. 
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S=S,-R, = 0.761-120° = 10,958mm’. 


2) 


Mathcad calculations give the following values: S=10,820 mm’, p,,,i, =9.167 mm. As can be seen in this 


example, empirical formulas accurately predict the basic parameters of the geometry of the end section of the gerotor 
machine. 

Discussion and Conclusion. The conducted research was aimed at solving the problem of improving the working 
bodies of gerotor hydraulic machines. It sufficiently corresponds to the works of domestic and foreign authors [10-15]. 
The key difference is that when constructing the rotor and stator profiles, ideal hypocycloids are used, providing the 
greatest kinematic accuracy of movement and contact interaction of the engagement elements. It was proposed to solve 
the well-known problem of providing tension in the coupling with respect to the difference in the parameters of the 
roller radius. A special algorithm was developed to control the smoothness of the profiles [12]. The optimal geometry of 
the profile, as in [10], was proposed to be determined on the basis of the maximum open area of the working body of 
the gerotor machine. But at the same time, the reduced contact radius should also have the largest value, as a guarantee 
of the lowest contact stresses. In this study, based on a large number of calculations using a purpose-built Mathcad 
program, in the range of rotor teeth z, from 2 to 7 and contour diameter R, from 100 to 200 mm, optimal values of the 
equidistant hypocycloid distance, or roller radii R., at which the best combination of the reduced contact radius and the 
open area of the working body of the gerotor was achieved. It turned out that the introduction of relative values of the 
radius of the roller, the open area, and the reduced contact radius provided obtaining absolute values of the same 
magnitudes for any contour radius within the established ranges using simple empirical dependences at an early stage of 
design. 

The developed technique makes it possible to evaluate the performance and operability of the gerotor machine at the 
design stage of the rotor and stator profiles. In the future, similar studies can be carried out for epicycloid and mixed 
engagements in gerotor hydraulic machines. 

The research results can be used in mechanical engineering under the design of gerotor hydraulic machines to 
improve their operating characteristics. 
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Abstract 
Introduction. Long-stroke movements in automated pneumatic drives account for a significant number of executive 


movements in coordinate tables, automated warehouses, cutting machines, etc. Long-stroke movements degrade the 
dynamic quality and positioning of the drive. This is due to the friction of the piston and the nonlinear characteristics of 
the compressed gas flow in significant volumes of the pressure and drain cavities of the cylinder. Thus, it seems 
promising to create an automated position pneumatic actuator for long-stroke movements. This will increase the 
productivity of processes while providing the declared accuracy. The objective of the work is to obtain a mathematical 
model and dependences of the critical parameters of the proposed automated position long-stroke pneumatic drive of 
fabrication system in the areas of acceleration, steady-speed movement, deceleration, and braking. 

Materials and Methods. The basis for calculations and modeling was the scheme of two trajectories of movement from 
point A to point E, taking into account the forces expended on these processes. The optimal displacement was 
determined using the Portnyagin’s principle (i.e., optimal performance). Proportional drive control was presented as a 
method of achieving the result. For long-stroke drive movements, schematic solution and design scheme were 
visualized in detail (presented as drawings). An original jet sensor with an internal pneumatic connection and a 
pneumo-mechanic discrete-proportional device for the control loop performance were proposed. The mathematical 
model included the movement and braking of the piston, the balance of mass flow, the pressure at points, and the 
control loop. The system of equations was solved by the Runge-Kutta method in the SimInTech software product. 
Based on the results of the study of a generalized mathematical model, the dependences of changes in the kinematic, 
power and pneumatic properties of the drive were constructed in real time during a typical positioning cycle. The 
information was summarized and presented as a set of graphs. 

Results. The mathematical model was formed according to a set of calculations. It took into account the dependences 
characteristic of the movement of the piston of the pneumatic cylinder. The balance of mass flow was investigated by 
the equations of gas flow during compression in the chamber, through distributors and chokes, in the discharge and 
drain cavities and in the control device. Inequalities describing the pressures at the points and the control loop were 
considered. A complex mathematical model was solved in the SimInTech software environment by the Runge-Kutta 
method with a variable integration step. A fragment of the program was selected as one of the illustrations. It showed 
that the software used the following indicators for calculations: target and reduced coordinates; absolute gas constant; 
coefficients of spring stiffness, resistance, adiabatic and viscous friction in the piston; compressor pressure; mass of the 
moving parts of the pneumatic actuator; strength of external resistances; diameters of the pipeline, the pneumatic 
cylinder piston and the braking device; length of the stroke of the cylinder piston; area of piston cavities and throttles; 
length of the pipeline and its internal volume. Thus, the program manipulated a significant set of data, which made it 


possible to obtain meaningful and adequate results. The relationship of blocks and diagrams used in solving the model 
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was schematically shown. We are talking about graphs of movements, areas, pressures, velocities and temperatures. 
Blocks with the program text and intended for integration were used. Thus, a mathematical model of an automated 
pneumatic drive of the fabrication system and the dependences of the basic parameters of its operation were obtained. 
The graphs indicated that the operating mechanism of the pneumatic actuator properly followed the proposed trajectory. 
Discussion and Conclusion. The research results allowed us to consider several stages of long-stroke movement of the 
drive, to determine the time frame of these processes (from 0 to 0.65s), as well as changes in pressure and speed of 
movement of the pneumatic cylinder carriage recorded in these intervals. There were five such stages: acceleration, 
steady-speed movement, deceleration, movement with positioning speed, and braking. Further research will focus on 


optimizing the system to reduce the duration and maintain accurate positioning under external influences. 


Keywords: long-stroke pneumatic actuator, jet control system, mechatronic module, pneumatic sensor, pneumatic 


actuator positioning, SimInTech software environment, Runge-Kutta method 
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AHHOTalna 

Beedenue. Ha mHHOXOOBEIe NepeMeLIeCHHA B ABTOMATH3MPOBaHHbIX THEBMONPHBOaX MPHXOAMTCA 3HAYMTEIIbHOe 
KOJIM4eECTBO HCHOJIHUTCIbHBIX JBHKCHHM B KOOPAMHAaTHbIX CTOIAaX, Ha ABTOMATH3HPOBaHHbIX CKIayax, pacCKpOHMHBIX 
MalmMHax UT. 1. J MHHOXOJOBbIC NepeMeleHHA YXyMMaloT WHHAaMMYeCKOe KayeCTBO H NO3HIMOHUpOBaHHe MpuBora. 
3ro obycioBieHo TpeHHeM MOpWHA HU HeMHeMHbIMM XapaKTepHCTHKaMH TOTOKa CxKaTOrO ra3a B 3HAYMTeJIbHBIX 
oObeMax HanopHoli U CIMBHOM MoNocTeHi WuMH Apa. Takum oOpa30M, ipecTaBsiaeTCA MepCHeKTHBHBIM CO3jaHve 
AaBTOMATH3HPOBaHHOTO NO3HIMOHHOTO MHEBMOMPHBOAa WIA JVIMHHOXOAOBBbIX TepeMelleHHH. ITO MO3BOJIMT TOBbICHTb 
TIPOW3BOMTeIbHOCTb TpoyeccoB mpu obOecneyeHHH 3aaBIeHHOH TouHocTH. Ilemb padoTbl — nomy4enne 
MaTeMaTHYeCcKOH MOeIM U 3aBHCHMOCTeli OCHOBHBIX THapaMeTpOB [peIO%KeHHOTO aBTOMATH3HpOBaHHOrO 
MO3HUMOHHOFO JVIMHHOXOAOBOrO MHEBMOMpHBOa TeEXHOOrMYecKorO OOOpyAOBaHHA Ha y4YacTKax pa3roHa, IBWKCHHA 
C YCTAHOBUBINeHCA CKOPOCTbH, 3aMeJICHHA UW TOPMOXKeEHHA. 

Mamepuanoi u memoooi. ba30ii 1a pacueToR HM MOJeMpoOBaHMA CTayia CXeMa J[ByX TpaeKTOpHi MepeMelleHHA V3 
touxw A B TouKy E c yyeToM CHI, 3aTpadeHHbIX Ha 9TH Mporecchl. OntTHMabHoe MepeMellleHve onmpeyemMM c 
TMOMOWbIO upHHuwna IloptasaruHa (TO eCTb ONTHMasIbHoro OpbicTpoyelictBua). I[ponopuMoHasbHoe yiupaBleHve 
IIPHBOJOM TipesCTaBsIeHO KaK MeTO AOCTWKeHHA pesybTaTa. J[y11 WMHHOXOAOBBIX NepeMelleHui MpHBosa TeTabHO 
BH3YasIM3HpOBaHbI (MpeCTaBIeHbl KaK PMCYHKH): CxXeMaTMYecKOe pellieHHe uM pacueTHad cxema. IIpenox%KeHbI 
OPHTHHAIbHbIM CTpyiHbIi WaTaMK C BHYTpeHHeli THEBMaTHYUeCKON CBA3bIO HU THeEBMOMeXaHW4uecKoe JMCKpeTHO- 
TIpOMOpWHOHAasIbHOe YCTpONCcTBO Wa ObicTpoyelicTBHaA KOHTypa yupaByeHua. MatemaTuueckasd MOeuIb BKIOUaeT 
TBYWKCHHe HU TOPMOXKeHHE MOPUWIHA, OamaHc MaCCOBBIX pacxOOB, TaBIeHHe B TOUKaX HM KOHTyp yupaBleHua. Cuctemy 
ypaBHeHu pelwanu Metoqom Pyure — Kytrst B uporpamMMHOM mpogzyKte «Cumuutex» (Simintech). Ilo utoram 


YMCCICTOBAaHuA odo6ujeHHOH MaTemMaTHueckon MOJCIH WOCTpOWJIH 3aBHCHMOCTH HU3MCHCHHA KHHCMaATHUCCKHX, 
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CHJIOBBIX HM IIHe€BMaTHYeCCKHX CBOMCTB TIpHBOJa B peaIbHOM BPeCMeHH IIpH THMOBOM UMKe MO3HIMOHHpOBaHHaA. 
Vndopmauto CyMMUpoBasIN HM IpeCTaBMJIM Kak COBOKYHHOCTb rpadHKos. 

Pe3yismamoi ucciedoeanua. Matematuueckas MOjesIb ctbopMupoBaHa MO KOMIWIeKcy pacueToB. Ona y4HTbIBaeT 
3aBHCHMOCTH, XapaKkTepHble JIA TBWKeHHA NOpWwHA WHeEBMOUMIMHApa. banaHc MaccoBbIX pacxOJOB MccyeqyerTca 10 
ypaBHeHHAM pacxojla ra3a Ip CxKAaTHH B KaMepe, Yepe3 pacipeyeuuTeIM U ApocceNu, B HarHeTaTeIbHOM UM CIMBHOM 
MOJIOCTAX HM B YIIpaBJIAIOWeM ycTpolicTBe. PaccMOTpeHbI HepaBeHCTBA, OMMCbIBaIOLIMe JaBIeHHA B TOUKAX MH KOHTYp 
ylpapienua. ClomwkHad MaTeMaTH4ecKad MOJeJIb pelllaacb B IporpaMMHoii cpeye «CumMuHTex» (Simintech) Metoyqom 
Pyure — KytTpi c u3MeHAeMBIM WaroM MHTerpupoBaHusa. DparMenT paooTsl NporpaMMbl BbIOpaH B KayecTBe OJHOL 
3 WWNOCTpalHi. OH MoKa3bIBaeT, 4TO COT 3afelicTByeT WIA pacueTOB TakKHe [OKa3aTeJIM, KaK: 3ajaHHad 
TIPHBeeHHbIe + KOOPAMHATbI; YHHBepcaylbHadt Ta30Bad MOCTOAHHaA; KOIPMUIMeHTHI %KeECTKOCTH Upy?KHHbI, 
COMpoOTHBIeHHA, ayuaOaTbI WH BA3KOTO TpeHHA B MOpWHe; WaBsleHwe KOMIIpeccopa; Maccy HOABIDKHbIX uacTeli 
MHEBMOMpPHBOa; CHJIy BHCWIHUX COMNPOTHBJICHHH; JHaMeTpbI TpyOonpoBoya, NOPWIHA MHEBMOLWMIMHApa HM TOpMO3HOTO 
yCTpolicTBa; MpOTAKeCHHOCTb XOJa NOPWHA WMIMHApa; MIOWaqWH MOpwHeBbIxX WoMocTeHi WH Apoccenel; WIMHy 
TpyOompoBoya MU ero BHYTpeHHHi oObeM. TakuM oOpa30M, MporpaMMa OlepupyeT 3Ha4HTeJIbHbIM KOMIIJICKCOM 
JAHHBIX, 4TO aeT BOS3MOXKHOCTS MOJYANTb CyU[CCTBCHHbIe HM aj[eKBaTHble pe3ybTaTbl. CxeMaTHYecKH MoKa3aHa 
B3aMMOCBA3b OJIOKOB HM WHarpaMM, HCIIOb30BaHHbIX IPH peleHHH Moyen. Peub ugeT o rpacdbuKax mepememjenui, 
TWIoWaye, TaBweHuli, ckopoctei u Temnepatyp. Ucnomp30BaHbl OOKH C TeEKCTOM IIporpaMMBI HM peqHasHayeHHble 
JUId WHTerpHpoBanna. TakuM oOpa30M MOJy4eHbI MaTeMaTHYeCKad MOJesIb ABTOMATH3MpOBaHHOTO MHeEBMOIpPHBOLa 
TeXHOJOrMYecKorO OOOpyAOBaHHA HU 3ABHCHMOCTH OCHOBHBIX WapaMeTpos ero paootur. [papukKn cBUeTeIBCTBYIOT O 
TOM, 4YTO HCHOMHUTCIBHbIM MCXaHH3M WHEBMONPHBOLa JOJDKHbBIM OOpa30M CylelyeT IpeA102%KCHHOM TpaeKTOpHH. 
O6cyatcdenue u 3akmiouenue. Vitorm padoTbl M03B0JIAIOT PaCCMOTPeTb HECKOJbKO 9TANOB [JIMHHOXOZOBOrO 
TIe¢peMeLIeHHA MIPHBOAa, ONpeAeIUTb BPeEMeHHbIe paMKU 3THX Uporeccos (oT 0 Ao 0,65 cek), a TakxKe (UKCHpyeMBIe B 
JaHHbIC IPOM@XKYTKH H3MCHEHHA aBICHHA HM CKOPOCTH ABWOKeHHA KapeTKH MHeBMOWMIMHApa. Takux 9TaloB MATS: 
pa3roH, JBYWKeHHe C yCTaHOBMBIIelCA CKOPOCTbIO, 3aMeJWICHHe, ABYXKeHHe CO CKOPOCTbIO NO3HUMOHMpoBaHHA 
TopMoxenue. JlaibHeliuMe vccieqoBaHua OyAyT cocpeyoOTOueHbI Ha OMTHMMH3alMH CHCTeMbI JIA coKpallleHHa 


TIPOAOJDKUTCJIBHOCTH HU NOATepKaHiaA TOUHOTO MOSHUMOHHPOBAaHHA TIpH BHCHIHHX BO3J]CHCTBHAX. 


Koroyesbie = CuI0Ba: ~=VIMHHOXOOBOH MHeBMONpHBOT, cCTpyiiHad cCHcTeMa ylpaBJIeHHs, MeXaTpOHHBbIM MOJYJb, 


TIHeEBMaTH4ecKHih AaATANK, MOSHUMOHMpPOBaHHe MHEBMOTIPHBOIa, IpOrpaMMHaad Cpeya «CHMHHTEX>», MeTO], Pyure — KyTrbi 


Baarojapnoctu: aBTOPbI BbIpaxKaroT OlaroqapHOcTE petakWHu UW peleH3eHTaM 3a BHHMAaTeCJIBHOe OTHOMICHHE K CTaTbe 


VW 3aMeudaHHA, KOTOPbIe MO3BOJINJIN MOBbICHTb Ce KadeCTBO. 


Aaa qWuTuposBanua. Kopotsr JI.A., Cunopeuko B.C., puxogzpKxo CI. Uccneqospanue WuHamuueckux XapakTepHcTHK 
aBTOMATH3UPOBaHHOTO NO3HIMOHHOFO JJIMHHOXOAOBOrO HHEBMONpHBOAa TeXHoOrMyecKorO OOopyzoBaHHA. Advanced 
Engineering Research (Rostov-on-Don). 2023;23(3):283—295. https://doi.org/10.23947/2687-1653-2023-23-3-283-295 


Introduction. The performance of the drives is determined by the accuracy of positioning and the speed of 
movement of the coordinates in different operating cycles. Modern fabrication systems are often equipped with 
automatic pneumatic actuators, which are characterized by long-stroke movements. These are, e.g., gantry resistance 
welding machines, coordinate tables, and cutting devices. 

Modern positioning pneumatic actuators for long-stroke movements in fabrication equipment provide a speed of 
up to 30 mm/s and an accuracy of ~1 % of the travel length. Customized actuators provide positioning accuracy of 
0.4 % at speeds up to 100 mm/s. Note that the trajectory of movements is formed by the compressed air flow control 
in pressure or drain pipelines and the pneumatic cylinder sides. In long-stroke drives, the length of such sides reaches 
3 m. Complex thermodynamic processes and compressibility in air flows are the key factors limiting the increase in 
accuracy [1-3]. 

Thus, it is necessary to increase the productivity of the working and engineering processes of the equipment while 


providing the declared accuracy. In this case, it seems promising to create an automated positioning pneumatic actuator 
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for long-stroke movements. The new solution should take into account such characteristics of the pneumatic actuator as 
speed, weight-size parameters, explosionproof, and fire protectability [2, 4]. 

The objective of the work is to obtain a mathematical model and dependences of the key parameters of the proposed 
automated positioning long-stroke pneumatic drive of industrial equipment in the areas of acceleration, steady-speed 
movement, deceleration, and braking. 

Materials and Methods. Figure | shows schematically the transport problem of moving from point A to point E 
along two trajectories. The forces expended on each of the movements are taken into account. The suboptimal 


movement ABCDE (trapezoid) with a simple control algorithm is realized in time 7, >> 7, > min. 


The optimal displacement AFE (bell) was obtained by solving the optimal speed based on the Portnyagin’s principle 


T,, > min, AL < |AL nox | . The result was achieved with more complex proportional control of the drive. The trajectory 


rt 


of movement provided the accuracy of switching motion controls along path L,,. 


Fig. 1. Movement trajectories: 1 — suboptimal ABCDE; 2 — optimal AFE 


In the description of trajectory 1, the switching points are indicated in Latin letters: A — for acceleration of the 


drive; B — for deceleration; C — for positioning speed; D — for stopping. In sections AB and BC, initial acceleration 


and braking are provided up to the positioning speed V,, and further stopping by the braking device AL, < A ae : 
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When moving along the second trajectory, A — switching to acceleration of the drive; F — switching to stop. 

An original jet sensor with an internal pneumatic connection and a pneumo-mechanical discrete-proportional device 
are proposed, which provides increasing the control circuit speed, since feedback in known analogues reduces the 
accuracy of the main engine by about 10-15 % during long strokes [4-6]. 

A schematic solution of a pneumatic positioning drive for long-stroke movements is shown in Figure 2. The drive 
operates in accordance with the suboptimal motion trajectory determined in the optimal speed problem for a given 
positioning accuracy. Here, IIL]1 — rodless pneumatic cylinder of long-stroke movements, which carries out the main 
motion; IIL2, 1113 — brake pneumatic cylinders that fix the drive; CA — jet sensor that determines the coordinate of 
movement, acceleration of the drive, its speed and force; P1 — pneumatic distributor with electropneumatic control, it 
controls the supply to the jet sensor; P2 — main control distributor; P3 — distributor with pneumatic control, it controls 
the operation of pneumatic brake cylinders; [1-4 — mufflers responsible for pressure relief into the atmosphere; 
Jl] — pressure sensor receiving data from the jet sensor; IIJIK — logic controller; LJ] — stepper motor controlling 
the distributor spool; bIIB — air-preparation unit; JJ[P1, J[P2 — throttle with check valve, used to regulate the speed of 


the rodless pneumatic cylinder of long-stroke movements of the main motion [7]. 
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Fig. 2. Circuit design proposal of an automated positioning long-stroke pneumatic actuator 


The drive contains a control system that monitors the position of the carriage of a rodless pneumatic cylinder, slows 
down when approaching the specified coordinates, and sends a signal to the braking gear [1, 4, 7]. 

The speed is increased by the introduction of a discrete proportional control device. Signals generated by the control 
loop are used for control. The device is made in the form of two nozzles and contains compensation measurements. 

Research Results. The design scheme of the pneumatic drive of long-stroke movements is shown in Figure 3. 

Here, IIL{1 — rodless pneumatic cylinder of long-stroke coordinate movements, carrying out the main motion; IL2, 
IIL3 — brake pneumatic cylinders that fix the drive during a stop in the desired position; CA — jet sensor that 
determines the coordinate of movement, acceleration of the drive, its speed and force; P1 — pneumatic distributor with 
electropneumatic control, it controls the supply to the jet sensor; P2 — main control distributor; P3 — distributor with 
pneumatic control, it controls the operation of brake pneumatic cylinders; [1-I'4 — mufflers responsible for pressure 


relief into the atmosphere; J[P1—/[P2 — throttle with a check valve, which serves to regulate the speed of a rodless 
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pneumatic cylinder of long-stroke movements of the main motion; JJ] — pressure sensor receiving data from a jet 
sensor; S — piston area of a rodless pneumatic cylinder of long—stroke coordinate displacements; P1—PS — 
investigated pressures at points 1-5; T1-T5 — investigated temperatures at points 1-5; F:» — friction force in a rodless 
pneumatic cylinder of long-stroke displacements; F,, — viscous friction force in a rodless pneumatic cylinder of long- 
stroke displacements; Fs, — external force in the rodless pneumatic cylinder of long-stroke movements; x — travel of 
the carriage of the rodless pneumatic cylinder of long-stroke movements; V — travel speed of the rodless pneumatic 
cylinder carriage; Crp — spring rate of brake pneumatic cylinders; m — mass to be moved; P; — pressure in brake 
cylinders; Py — pressure in the control channel; f1—/4 — areas of the through sections; Pa — atmospheric pressure; 
d3\1—d33 — diameters of distributor spools; Cappi—Cnpp3 — spring rate of distributors 1-3; xpi—xp4 — displacement of 
distributor spools 1-4; V;1-Vp4 — travel speed of distributor spools 1-4; Fyyi—F's42 — force of the distributor control 


electromagnet 1—2; P,, — compressor pressure; 7, — temperature generated by the compressor. 


Ms p4 


Fig. 3. Diagram of the automated pneumatic drive of long-stroke movements of the industrial equipment 
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The mathematical model was formed with the following assumptions [8-12]: 

— pressure of the compressed air source remained constant over time; 

— thermodynamic process of gas behavior in a pneumatic system was adiabatic; 

— description of pneumatic devices used the ideal gas model, since the pressure of the pneumatic system did not exceed 10 bar; 

— leaks were not taken into account; 

— viscous friction force was proportional to the velocity; 

— expense ratio was recognized experimentally through identification; 

— mass of the moving part was constant; 

— force at the output link of the pneumatic drive was constant. 

1. Equation of motion of the pneumatic cylinder piston [1, 4]: 

=S4(p)-7.)-Fy kay Fug sign a: F, (1) 
Here, S — area of the piston of the discharge and drain sides of the rodless pneumatic cylinder of long-stroke coordinate 
movements of the main motion, m?; P, P, — air pressure in the discharge and drain sides of the pneumatic cylinder, Pa; 


F, 


"ay —— external forces, N; ks: — viscosity friction coefficient; i — friction force, N; o — Boolean parameter: a = 0 at 
Ds < Pon, Ad a= 1 at p, > p,,.,; ps — pressure in the control channel, Pa; p,,— atmospheric pressure, Pa [12-14]; 
m — mass of the moving parts of the drive, kg; Fr — braking force, N. 
Ee a pH . Crp m* (255 ~ Xn ) ? (2) 
where Cnp + — spring rate of the pneumatic cylinder of the brake; ,1 — friction ratio. 
2. Equation of motion of the piston of the braking pneumatic cylinder: 
d?x dx 
m = = Cy m : x m Xin Sin m 1 ee m Kam m : 7 * (3) 
di” (% ) fi dt 
Here, Cup r — spring rate of the brake pneumatic cylinder; xo, — coordinate of the initial compression; S; — effective 
area of the piston of the drain side of the brake pneumatic cylinder, m?; P— air pressure, respectively, in the discharge 


side of the brake pneumatic cylinder, Pa; F's 7— external forces, N; ky: — viscosity friction coefficient. 
3. Mass expenditure balance equations: 


Ga (t)— Gay (1) + Gane (t) = 9: (4) 

G.,,. (t)—G,; (t)-G,, (t) + G.,, (t) =0, (5) 

Gays (t)— Gop (4) — Gyr (1) + Gage (t) = 0, (6) 

Giys (t)— Gaps (4) — Gaya (t) + Gag (1) = 0, (7) 

Cry (t)— Gps (4) + Ge (1) = 0 - (8) 

Here, G,,(¢) — gas mass flow rate under compression in the chamber; G,,(t),G,,(¢)— mass flow through 


distributors; G 


mi (t) and G,,,.(¢) mass flow rate in the discharge and drain chambers of a rodless pneumatic cylinder of 


long-stroke movements; G,,(t), G,,,(t), G,,.(t), G,,3(t), G,,(t) — mass flow through throttles in the drain line, at the 
inlet to the nozzle unit of the jet sensor and at the outlet of the nozzle unit; G,,,(¢), G,,.(4), G,, (0, G,,.4(4) — mass flow 


rate in the control channels of the control device, at the outlet of the control device, distributor of brake pneumatic 
cylinders [7]. 
awViod 
(=P ) 
4-E-R-T dt’ 
where p= 1.288 kg/m? — air density; V — volume of the chamber; R = 287 J/(kg-k) — gas constant; E — volume 


modulus of elasticity of air; T — temperature at the point; “ — pressure variation at the point [6]. 
t 
W, dp 
G,,, (t)=— - + 10 
mi (1) = k-R-T dt’ 
W, dp 
G,,,. (t) =—+—: =, 11 
mo (t) k-R-T dt oo 
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where W, and W, — current volumes in the pressure and drain sides of the base pneumatic cylinder, m+; & — adiabatic 
exponent (for air, k= 1.4). 
4. Pressure equations in points: 


ap, BART [p,, oe {2 va Top? PL a 12 
dt S(x+Xq, )JE, (x+X, ae ee ae 


k-1 


k- f.-.JR-T. 2k F 
Ps fa BT = Se ae (13) 


ie. dt” 


dt S(L-x+Xy) JE, P, 

dp, hv EST (2) ra 
= W “VP TE a ae rae (14) 
dp, _ k- | R T, k- Ss: R T, Ps —— 

PE MEE 


kK: R-T, 
dps _ ff = aA 


dt W,- 


Here, k — adiabatic exponent; R — gas constant, J/ — ; Tw — air temperature in the main line, K; Pa — 


atmospheric pressure, Pa; P; — Ps — pressure in the flow parts of the pipeline, Pa; W; — Ws — volumes of the flow 
parts, m?; €;— &7 — resistance values in the line; fi —_fs— areas of the flow sections of the pipeline, m’; L — maximum 
stroke of the piston, m; xo1, x02 — ratio of the initial volumes of the pneumatic actuator to the useful area of the piston 


side of the pneumatic cylinder, m; — piston speed, m/s. 
t 


5. Control loop equations. 
Equation of motion of the distributor spool 1: 


M, » i Cpa ta Fa sign( Fi (2 Fa (17) 

Equation of motion of the distributor spool 2: 
pp = eg Sn —P -sign( Se Fs (= Fyz: (18) 

Equation of motion of the distributor spool 3: 
ea 

Equation of motion of the distributor spool 4: 


Here, S,,— area of the distributor spool end face, m*; P3, Ps, Py — pressure in the control channels, Pa; F’.— resistance 
forces to the movement of the distributor spool, N; Fy — reaction forces of stops, N; F'sm — force of the electromagnet 
acting on the distributor spool, N; cup p — spring compression ratio, N/m; m;, — weight of the distributor spool, kg. 

To solve the mathematical model, SimInTech software product was used. We applied the Runge—Kutta method with 


a variable integration step (Fig. 4—5). 
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Fig. 4. Model fragment in SimInTech program 


Here, P,1—P,3 — blocks with the program text; f1, f2 — graphs of the output of areas | and 2 of the jet sensor 
throttle; fis: — block of the problem of the areas of the flow sections of the jet sensor; fis2 — graphs of the output of the 
areas of the flow sections of the jet sensor; d71-d7T8 — integration blocks; P1—P4 — graphs of the output of the 
received pressures at points 1-4; P — general schedule for the output of all pressures; V — graph of the obtained speed 
output of the rodless pneumatic cylinder carriage; x — graph of the received movement output of the rodless pneumatic 
cylinder carriage; x, V — general graph of the output of the received movement and the speed of the rodless pneumatic 
cylinder carriage; 7 — graphs of the output of the obtained temperatures; ¥3 — graph of the received movement output 
of the control distributor spool. 

2 Se Re nperpagoo an L=0.43 
ee S=(pi*D*2)/4; 


QA *eOx|1|Se/At|E FUR elrv : 
1 input pl, p2, p3, p4, V, x,fs1, fs2, ds; St=(pi*Dt*2)/45 


output pl, p2, p3, p4, V, x, d, T1, 12, T3, T4, TS, T6,p5; | kt=0.4 
init p1=100000,p2=100000, p3=100000, p4=100000, V=0, x=0; begin 
xk=0.31 30 | if x<xk+@.0577 then 

- R=287; 16=0 
Tm=293; foe 
ks1=653 else 
ks2=100; p6=pm 
ks3=70; end 

10 ks4=70; 5 begin 
d3=0.001; ; . 
d5-0.006; if x>=0.363 then 
pm=6*1045; p5=pm 
pa=1.01325*10°5; else 
k=1.45 p5=0 
x01=0.008 ae | end 
02=0.008; 

Acs N=p5*St; 
=6 5 
Fvn=20; Ft=kt*N; 

20 dtr1=0.005; begin 
dtr2=0.005; if x<xk then 
peel d2=0.008 
Dt=0.04; 1 
kvt=320; bag 

d2=0.0004 
end; 


b) 
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“begin a 
: begin 
so | if x<xk+@.0577 then 4f p3>=600000 then p3=600000; 
d=ds = 
1 5 begin 
i 7 a lls then p3=pa; 
— = en 
s 20 begin 
end; if p4<=pa then p4=pa; 
© begin end 
if x<xk then T= (p/p) *( (K-27) *Tm5 
5-0 T2=(p2/pm)*((k-1)/k)*Tm5 
a T3=(pm/pm)*((k-1)/k)*Tm; 
else La ila a 
i T5=(pm/pm)*((k-1)/k)*Tm; 
5=(pi*d5*2)/4 T6=(p4/pm)*((k-1) /k)*Tm5 
60 end; G11=(k*f1*sqrt(R*Tm) /(S*(x+x01)*sqrt(ks1) ))*sqrt(abs(pm*2-p1*2) ); 
1=(pi*d*2)/4; 20 G12=(k*f3*sqrt (R*T3) /(S*(x+x@1)*sqrt(ks2)))*(p3/p1)*((k-1)/(2*k))*sqrt(abs(p1*2-p3*2)); 
- - G13=(k*f3*sqrt(R*T5) /(S*(x+x01)*sqrt(ks2) ))*(p4/p1)*((k-1)/(2*k))*sqrt(abs(p1*2-p4*2)); 
#2=(pi*d2*2)/4; G14=(k*p1/(x+x01))*V; 
£3=(pi*d3*2)/4; G21=(k*f2*sqrt (R*T2) /(S*(L-x+x02)*sqrt(ks2) ))*(p2/pa)*((k-1)/(2*k))*sqrt(abs(p2*2-pa%2)); 
G22=(k*p2/(L-x+x@2) )*V; 
Ltr=0.01; G31=(k*#3*sqrt (R*T3) /(W1*sqrt(ks3)))*sqrt (p1*2-p342); 
cpr=500000; G32=(k*fs1*sqrt(R*T4) /(W2*sqrt(ks4)))*((p3/pa)*((k-1)/(2*k)))*sqrt(abs(p3*2-pa*2)); 
ee eee G41=(k*f3*sqrt (R*T5) /(W1*sqrt(ks3)))*sqrt(p1*2-p4%2) ; 
W2=Ltr*((pi*dtr2*2)/4); G42=(k*fs2*sqrt(R*T6) /(W2*sqrt(ks4)))*((p4/pa)*((k-1)/(2*k)))*sqrt(abs(p4*2-pa*2)) 5 
iS begin p1'=G11-G12-G13-G14; 
ey 100 p2"=-G21+G622; 
if L<x then Fy=(cpr*(x-L)) p3"=G31-G32; 
70 | else p4'=G41-G42; 
if (L>=x) or (x>=0) then Fy=0 ae ‘io ogeticy haciad epee 
else if x<@ then Fy=cpr*x SR a ea rr mee, fe nea eye rahe 
c) d) 


Fig. 5. Programming block in SimInTech: a — part of the source data; b — part of the logical functions; 


c — part of the assignment of variables; d — skeletal code 


The study of the generalized mathematical model of the proposed drive made it possible to obtain graphs of the 
behavior of the drive during acceleration, deceleration, and stop (Fig. 6), describing the changes in kinematic, power 


and pneumatic properties of the drive during a typical positioning cycle in real time [15]. 


P5, Pa V,m\s x,m P4,Pa P3,Pa P2,Pa P1,Pa V,m/s P3,Pa P4,Pa x,m P2,Pa Pl, Pa P5, Pa 
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Fig. 6. A set of graphs based on the data of a generalized mathematical model 


Discussion and Conclusion. The graph shows a long-stroke movement along the trajectory proposed in Figure 1. 
The operation of the drive consists of several stages. 

1. 0-0.1 sec. Acceleration. The pressure in the pressure side increases to 4 bar, the pneumatic cylinder carriage 
speed — | m/s. 

2. 0.10.38 sec. Movement at a steady speed. Pressure in the discharge side — about 3.8 bar. Pressure in the drain 


line increases to 1.6 bar, the pneumatic cylinder carriage speed — 0.85 m/s. 
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3. 0.38-0.5 sec. Slowing down. Pressure in the pressure and drain sides increases. The pneumatic cylinder carriage 
speed decreases to 0.075 m/s. 
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4. 0.5-0.65 sec. Movement with the speed of positioning. Pressure in the pressure and drain sides increases, the 
pneumatic cylinder carriage speed — 0.075 m/s. 

5. 0.65 sec. Switching to braking, activation of an external braking gear. 

The graphs obtained confirm that the long-stroke movements of the pneumatic actuator are performed in accordance 
with the proposed trajectory (Fig. 1), and the control system functions properly. Further research will focus on 


optimizing the system to reduce the duration and maintain accurate positioning under external influences. 
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Abstract 
Introduction. In recent decades, knowledge about DNA has been increasingly used to solve biological problems 


(calculations using DNA, long-term storage of information). Principally, we are talking about cases when it is 
required to select artificial nucleotide sequences. Special programs are used to create them. However, existing 
generators do not take into account the physicochemical properties of DNA and do not allow obtaining sequences 
with a pronounced “non-biological” structure. In fact, they generate sequences by distributing nucleotides 
randomly. The objective of this work is to create a generator of quasirandom sequences with a special nucleotide 
structure. It should take into account some physicochemical features of nucleotide structures, and it will be 
involved in storing non-biological information in DNA. 

Materials and Methods. A new GATCGGenerator software for generating quasirandom sequences of nucleotides was 
described. It was presented as SaaS (from “software as a service”), which provided its availability from various devices 
and platforms. The program generated sequences of a certain structure taking into account the guanine-cytosine (GC) 
composition and the content of dinucleotides. The performance of the new program algorithm was presented. The 
requirements for the generated nucleotide sequences were set using a chat in Telegram, the interaction with the user was 
clearly shown. The differences between the input parameters and the specific nucleotide structures obtained as a result 
of the program were determined and generalized. Also, the time costs of generating sequences for different input data 
were given in comparison. Short sequences differing in type, length, GC composition and dinucleotide content were 
studied. The tabular form shows how the input and output parameters are correlated in this case. 

Results. The developed software was compared to existing nucleotide sequence generators. It has been established that 
the generated sequences differ in structure from the known DNA sequences of living organisms, which means that they 
can be used as auxiliary or masking oligonucleotides suitable for molecular biological manipulations (e.g., amplification 
reactions), as well as for storing non-biological information (images, texts, etc.) in DNA molecules. The proposed 
solution makes it possible to form specific sequences from 20 to 5,000 nucleotides long with a given number of 
dinucleotides and without homopolymer fragments. More stringent generation conditions remove known limitations and 
provide the creation of quasirandom sequences of nucleotides according to specified input parameters. In addition to the 
number and length of sequences, it is possible to determine the GC composition, the content of dinucleotides, and the 
nature of the nucleic acid (DNA or RNA) in advance. Examples of short sequences differing in length, GC composition 
and dinucleotide content are given. The obtained 30-nucleotide sequences were tested. The absence of 100 % homology 
with known DNA sequences of living organisms was established. The maximum coincidence was observed for the 
generated sequences with a length of 25 nucleotides (similarity of about 80%). Thus, it has been proved that 
GATCGGenerator can generate non-biological nucleotide sequences with high efficiency. 


© Kiryanova OYu, Garafutdinov RR, Gubaydullin IM, Chemeris AV, 2023 


Kiryanova OYu, et al. GATCGGenerator: New Software for Generation of Quasirandom Nucleotide Sequences 


Discussion and Conclusion. The new generator provides the creation of nucleotide sequences in silico with a given GC 
composition. The solution makes it possible to exclude homopolymer fragments, which improves qualitatively the 
physicochemical stability of sequences. 


Keywords: GATCGGenerator, nucleotide sequences generator, synthetic nucleic acids, random sequences, data storage 
in DNA, steganography, NYRN-oligonucleotides, calculations with DNA, cryptography, DNA-tagging in hydrology 
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AHHOTauHA 

Beedenue. B nocneyuue AecaTuneTua 3HaHHa Oo JIHK sce wipe MpHMeHsAIOTCA JIA pelleHHA HeOMOMOTM4eCKHX 3aa4 
(BEIuuceHHA c MoMoubIO JIHK, JonroppemeHHoe xpaHeHHe uH@opMayHn). B nepBytlo O4epeyb perb UeT O CuyyasAx, 
Kora HeEOOXOAMMO TOAOOpaTb HCKyCCTBeHHbIe HYKICOTHAHbIe MOCIeAOBaTeIbHOCTH. JIA Ux CO3TaHUA MUCHOb3YIOTCA 
creljMasIbHble MporpaMMbl. OjHako cyljecTByrollMe reHepaTopbl He YAHTHIBaIOT (PU3HKO-xXuMHYecKHe cBolicTBa JIHK 
WM He MO3BONAIOT MOMWYUATb MOCIeAOBaTeIbHOCTH C ABHO BbIPAKEHHOM «HEOMOMOrHMYeCKOM» CTpyKTypoi. DakTH4ecKH 
OHH TeHEPpHpyIOT MOCIeAOBaTeNbHOCTH, pacnipeyeIAA HYKICOTHAbI CydaHbIM oOOpa30M. Llenb1o waHHoM padoTsl 
ABJIACTCA CO3aHHe TeHepaTopa KBa3HCIIyYaHHbIxX MOCIeAOBAaTeIbHOCTeli C OCOOOM HyKeOTHAHOM cTpyKTypou. Ou 
JOJDKCH YYHTLIBATh HEKOTOpbIe (PH3MKO-XUMHYeCKHe OCOOeCHHOCTH HYKJICOTHAHBIX CTpyKTyp MH OyseT 3alelicTBOBaH 
lip xpaHeHun HeOvoorMyeckol HHdopmauH B THK. 

Mamepuanot uu memodei. OrucaHo HoBoe porpaMMHoe oOecneyeHue GATCGGenerator ana reHepaluu 
KBa3HCIydalHbIx MOCIeOBaTeIbHOCTelHi HyKeCOTHAOB. OHO mpegoctaBiaeTca Kak SaaS (oT auru. software as a 
service — lIporpaMMHoe oGecrieyeHHe Kak ycylyra), 4TO OOecHe4MBaeT erO MOCTYMHOCTb C pa3HbIX yCTpOvcTB U 
tu1aTdpopM. IIporpamMa reHepupyeT NOcieqOBaTeIbHOCTH OMpeseIeHHOM CIpyKTypbl C yYeTOM ryaHHH-WMTO3HHOBOTO 
(GC) coctaBa u cogxepxKaHuaA WuHyKIeoTuAOB. IpeyzctaBeHa padota alropHTMa HOBO UporpamMpl. TpeOoBaHua K 
CreHepHpOBaHHbIM HYKJICOTHAHBIM TOCIICQOBaTeIbHOCTAM 3a/[aHbI C MOMOIIbIO YaTa B «Temerpam» (Telegram), 
HaraqHO 1OKa3aHO B3AaMMOJeiCTBHe C MOIb30BaTeIeM. OnpeyeseHbI H OOOOMeHbI pa3ssIM4HA BXOJHBIX MapaMeTpoB U 
TIOJY4aeMBIX B Pe3yIbTaTe paOoTh! NPOrpaMMBbI KOHKPeTHBIX HYKJICOTHTHBIX CTpyKTyp. TaloKe B COMOCTaBJICHHM JjaHbl 
BPeMeHHbIe 3aTpaTbl TeHepallMuH MOcIeAOBaTebHOCTeH Tp pa3IM4HbIX BXOJHBIX aHHbIX. V3yyenbr KOpoTKue 
TOCHeOBATEIbHOCTH, pa3ziMuarouMeca 0 Tuny, WuHe, GC-coctasy u coxepxKaHHIo AHHyKeoTHA0B. B Tabm“4HOM 
BH e TOKa3aHO, Kak B 9TOM Cilyuae COOTHOCATCA BXOJHBIC HM BbIXOJHbIe WapaMerTpsl. 

Pe3yivmamoi uccaedosanua. Co3yanHoe UporpaMMHoe obecneyeHHe cpaBHWJIM C CYLIeCTBYIOIMMU reHepaTopaMu 
HYKJICOTHIHBIX TMOCIe€LOBaTeIbHOCTeH. YcTaHoBsIeHO, 4YTO TeHepHpyeMble NMOCMeqOBaTebHOCTH OTIIMYAIOTCA 10 
CTpykType oT u3BecTHBIx JJ[HK-nocneqoBaTesbHOCcTeH %KMBbIX OPraHH3MOB, a 3HAYHT, MOryT ObITb HCHOJIb30BaHbI B 


KauyecTBe BCHOMOraTeCJIBHbIX HWJIM MaCKHpyrOUHxX OJIMTOHYKIJICOTHAOB, MPHrOWHbIx 1A MOJI€KyJIAPHO-OMOJIOTH4eCKHX 
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MaHHIyAui (HalpHMep — peakUHH aMiiMpuKkallun), a Take WJId XpaHeHua B MoeKysax JIHK HeOuonormueckon 
nudopMannu (H300paxKeHHH, TekcToBUT.y.). IpeqmoxeHHoe pelieHHe aeT BO3MO%KHOCTb (pOpMHpoBaTb 
cHelHpuyeckne MocteqOBaTebHOCTH AIHHON oT 20 yo 5 000 HykNeoTHAOB C 3aaHHbIM YHCIIOM JHHYKICOTHOB H 
6e3 rOMONOJIMMepHBIX y4acTKoB. bosiee *KeCTKHe YCIOBHA TeHepalluH CHUMAIOT H3BECTHbIe OrpaHHYeHHA UH MO3BONAIOT 
CO31aBaTb KBa3HCIydaliHble MOCICOBATeIbHOCTH HYKJICOTHOB TO 3aaHHbIM BXOJHBIM TapametpaM. Kpome 
KOJIMYECTBA H JVIMHbI NOCIeAOBATeIbHOCTeH MO2KHO 3apaHee OlpexzemuTb GC-coctas, coqepxKaHve JMHYKIICOTHOB U 
pHpoxy HykKieuHopon KucnoTer (HK usm PHK). I[pusogatca mpumepbt KOpoTKHX mMocseqoBaTesbHOCTeH, 
pa3zmuaromjuxca 0 WIMHe, GC-cocTaBy Mu coyepxKaHuIo AuHyKIeoTHA0B. Tlomyyenupie 30-HyKiIeoTHAHEIC 
MOCEAOBATeIbHOCTH NPOWIM MpoBepky. YcTaHoBeHO oTcyTcTBHe 100-npoyeHTHO romomoruu c w3BecTHEIMU JJHK- 
MOCJICHOBaTeIbBHOCTAMH %KMBbIX OpraHH3MoB. MakcuMasibHoe coBlayleHve HaOmroyaoch [Id creHepHpoBaHHBbIx 
TOceqoOBaTesbHOCTeH IMHO 25 HyKNeoTHyOB (cxogcTBo oKoO0 80%). Takum o0pa30M joKa3aHO, 4TO 
GATCGGenerator MoxeT Cc  BBICOKOM 9eKTHBHOCTbIO YTeHepHpoBaTb HeOMONOrM4YeCKHe HYKJICOTHAHBIC 
MIOCJI€OBATeIbHOCTH. 

O6cystcoenue u 3akiro4“enue. Hopbiit reHepaTop NO3BOAeT CO3TaBaTb HYKIICOTHAHbIe MOCIeOBaTeIbHOCTH in silico c 
3ayjaHHbim GC-coctaBom. PelmieHHe aeT BO3MO2%KHOCTb MCKMIIOUMTb TOMOMOJMMepHble (parMeHTbI, YTO KaYeCTBEHHO 


ystyaiaetT (pH3UKO-xXHMUyeCK yO cCTaOHJIBHOCTb MOCIIENOBaTeIbHOCTeH. 


Kiroueesie —cnoea: + GATCGGenerator, reHepaTop HYyKIIeCOTHZHbIX MOCIeAOBaTebHOCTeH, CHHTeTHUeCKHE 
HYKJICHHOBbIC KHCJIOTHI, CilyualiHble MOCIeTOBaTebHOCTH, XpaHeHHe JaHHEIxX B JIHK, creraHorpadua, NYRN- 


ONMPOHYKICOTHABI, BBIYHCIeHHA Cc MOMOobIO JIHK, kpuntorpadua, JIHK-merunkn B ruyqpomornu 
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Introduction. DNA is a unique biopolymer that provides storage, transmission and reproduction of genetic 
information in living organisms. DNA molecules consist of four types of nucleotides containing nitrogenous 
bases: adenine (A), guanine (G), cytosine (C), thymine (T). Their possible combinations provide nucleotide 
sequences forming functional genetic elements. In molecular biology and genetics, the basic investigations are 
carried out on nucleotide sequences of living organisms, but there is an increasing need to create artificial 
sequences, especially, when solving non-biological tasks (e.g., DNA calculations [1, 2], storage in DNA [3], 
cryptography [4], DNA tags in hydrology [5], etc.). 

It is expected that by the end of 2040, the volume of information will reach several yottabytes (1074), which requires 
its structuring and storage. Both of these processes affect significantly the consumption of energy resources, as well as 
the production of storage devices and peripheral devices (hard drives, solid-state drives). To store such an amount of 
information, more than 10° kg of extra pure silicon is required [6], which may not be enough. The solution is seen in 
using the principles of DNA to work with large-scale amounts of data. 

Nucleotide sequences are easily digitized by assigning the corresponding binary codes to individual nucleotides 
[7-11] or blocks of nucleotides [12—14]; therefore, text, graphic or multimedia files can be converted into nucleotide 
sequences [15-18]. Artificial nucleotide sequences can be made manually or generated using special software (DNA 


generators), depending on the tasks being solved. Some DNA generators were developed as independent applications, 


Kiryanova OYu, et al. GATCGGenerator: New Software for Generation of Quasirandom Nucleotide Sequences 


others — as part of software packages designed to solve general [19] '* 7+ 5 or specific tasks [20]. As a rule, DNA 
generators are developed on the basis of combinatorial approaches and produce random sequences of a given length of 
guanine-cytosine (GC) composition. However, such software solutions do not take into account the chemical properties 
of nucleotides and do not provide obtaining sequences with a certain structure (e.g., without homopolymer sites or long 
repeating motifs). Therefore, the sequences created by such generators cannot always be reproduced in the laboratory. 
Moreover, such sequences may be identical to DNA fragments existing in nature, which introduces ambiguity when 
trying to encode information of a non-biological nature. 

The presented work is aimed at creating a generator of nucleotide sequences of a special structure that can be used 
when encoding text, graphic and other information in DNA molecules. 

Materials and Methods. The criteria that should be kept in mind when creating sequences were defined. The need 
to vary the GC composition, set a certain number of dinucleotides, and exclude homopolymer sites in sequences was 
taken into account. 


A team of authors has developed the GATCGGenerator program in Python 3.6 (Anaconda distribution)°. To create a 
bot’ in Telegram, Numpy 1.19 [21] and the Python GATCGGenerator library were used. The solution was provided as 
SaaS (from “software as a service”), which opened up the possibility of access from different devices and platforms. 

Input parameters included the number of sequences, their length, GC composition, and dinucleotide content. The 
generator excluded repeats with a length of two nucleotides more than four times. The result was presented as a CSV 
file, which contained the following information: sequence, GC composition, and the number of all nucleotides. 

Repeats and homopolymer fragments were stored as a separate list. First, a sequence of four elements was randomly 
generated (random.choice(nuk), where nuc = 'ACGT'). Then the search for repetitions was performed. If there was at 
least one item from the list, a new random generation was performed. Next, the GC and NN composition was 
calculated. If the NN composition did not match the user-defined range, the paired nucleotide was replaced randomly 
and the GC composition was recalculated. If the sequence matched the input parameters, it was written to a set of 
sequences. 

Below is the operation of the program algorithm. 

Type, GCmin, GCmax — range of possible GC content, NNmin, NNmax — range of possible dinucleotide content 
NN%, N — quantity, S — sequence, 1 — sequence length, count — total number of sequences 
Pseudocode 
Start 
Input (Type, GC, NN, N) 

Comprehension of a list of repeating motifs, homopolymer sites rep.list 
Count = 0 
sequences = set() 
IFi<N? 
IF (rep.list(k) C S?) 
Return to step 1. 
ELSE 
NN = len(DL REGEX. findall("join(S))) 
NN_perc = (NN x 2/1) x 100 
IF NNmin <NN_perc < NNmax 
GC = S.count('G') + S.count('C') / 1 x100 
IF GCmin < GC < GCmax 
IF type == DNA 
Step 2. 
A_perc = S.count('A') /1 x 100 


' Nucleotide Sequence Generator. nucleotide-generator.herokuapp.com. URL: https://nucleotide-generator.herokuapp.com/ (accessed: 01.12.2022). 
2 DNA Sequence Tools: Random Sequence Generator. molbiotools.com. URL: http://www.molbiotools.com/randomsequencegenerator.html 
(accessed: 01.12.2022). 

> Random DNA Sequence Generator. faculty.ucr.edu. URL: http://www. faculty.ucr.edu/~mmaduro/random.htm (accessed: 02.12.2022). 

+ Random DNA Sequence GenScript. genscript.com. URL: https://www.genscript.com/sms2/random_dna.html (accessed: 04.12.2022). 

5 Random DNA Generator. Computer software. URL: http://54.235.254.95/cgi-bin/gd/gdRandDNA.cgi (accessed: 04.12.2022). 

® Anaconda. Anaconda Inc. anaconda.com. URL: https://www.anaconda.com/ (accessed: 20.01.2023). 


7 Python telegram bot. github.com. URL: https://github.com/python-telegram-bot/python-telegram-bot (accessed: 01.12.2022). 
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G_perc = S.count('G’) /1 x 100 
C_perc = S.count('C’) /1 x 100 
T_perc = S.count('T') /1 x 100 
U_perc = S.count('U') /1 x 100 
Count = count +1 
sequences.add(S) 
ELSE S = S.replace('T’, 'U') 
Step 2.) 
ELSE 
Return to step 1. 
ELSE 
Random replacement of the second repeated character, 
GC = S.count('G') + S.count('C') / 1 x100 
Output Sequences: (S, GC%, NN%, A%, G%, C%, T/U%) 
End 


The requirements for the generated nucleotide sequences were set using Telegram chat. An example of user 
interaction is shown in Figure 1. 


GATCGGenerator 
bot 


What can this bot do? 


Online generator of DNA/RNA sequences with specified GC and 
dinucleotides content 


August 26 


start 16:06 W 


Send a generator type j¢.9- 


909 


DNA 16:06 4 


DNA generator was chosen. Starting the process 4 ¢.p¢ 


( Enter the number of sequences to generate: 4 ¢.9¢ 
om 50 16:06 

") Enter the sequence length: 4 ¢.9¢ 

@ ‘50 16:06 w 


Specify the GC content in %. Use the following format: <min>- 
<max> 16:06 


50-60 16:06 W 


Specify the dinucleotides content in %. Use the following format: 
<min>-<max>. If dinucleotides are not mandatory send 0-0. 4 ¢.9¢ 


19-20 16:06 W 


Your job is queued. Please, wait for the result. 4 ¢.9¢ 


© 898 98 


Fig. 1. Example ofa user chat in Telegram 
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In the framework of the presented work, the functionality of random sequence generators and GATCGGenerator 
were compared. The differences between the input parameters and the specific nucleotide structures obtained as a result 
of the program were determined (Table 1). 


Table | 
Comparison of GATCGGenerator functionality to other nucleotide sequence generators 
DNA 
. Ri 
Nucleotide eeuuenee endon Random Random 
GATCGGenerator Tools: DNA 
Sequence DNA DNA 
[20] : Random Sequence i io 
Generator 19 | Sequence Generator 
Sequence Generator 
Generator? 
Maxi length 
ee 5,000 1,000,000 10,000 1,000 
(nucleotides) 
1; 10; 50; 
Number of sequences 100 1 100 100 
Input GC composition " é — + (*) 
(%) 
GC composition (%) number — number 
Input NN interval 
composition (%) — 
No homopolymer 2 
sites 
DNA/RNA / 
Sequence type DNA/RNA DNA ie : DNA 
Protein 
Output of results CSV file Text on the screen 
(*) User enters AT composition 


GATCGGenerator has a broader functionality, it allows the user to specify the number of dinucleotides, create 
sequences without extended homopolymer sites and repeats that affect the success of the experiment. In existing 
generators, it is only possible to vary the GC composition. 

The program created by the authors of this research generates a given number of quasi-random sequences of 
nucleotides that do not have homology with natural DNA, but are suitable for molecular biological manipulations. 

Research Results. GATCGGenerator allows you to generate specific DNA or RNA sequences from 20 to 5,000 
nucleotides long, containing a given number of dinucleotides and not containing homopolymer sites (no more than two 
identical nucleotides located side by side). More stringent generation conditions can cause a long selection of 
sequences. As an example, we give a small range of possible content of guanine and cytosine and dinucleotides (e.g., 
GC composition 45-50 % and NN composition 10—20 %). The operating period of the program for various input data is 
shown in Table 2. 


Table 2 
Sequence generation time for different inputs 
Data inputs . 
Time (s) 
Length Number GC, % NN, % 
20 10 50-60 20-50 3.45 
30 10 50-60 20-50 3.91 
20 10 50-60 40-50 9.74 
30 10 50-60 40-50 9.53 
30 10 40-50 20-20 8.80 
1,000 100 45—50 40-50 11.49 
2,000 100 45—50 10-20 240.25 
5,000 100 50-60 20-50 11.57 
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GATCGGenerator, through more stringent sequence generation conditions, removes the limitations of known DNA 
generators and creates quasi-random sequences of nucleotides depending on the indicated input parameters. You can 
specify the required number of sequences, their length, GC composition and dinucleotide content, as well as the nature 
of the nucleic acid (DNA or RNA). Specifically, sequences created by GATCGGenerator can be used in DNA 
steganography, applied to protect and transmit information through hiding the message content in the nucleotide 
sequence [3]. 

The proposed software solution (GATCGGenerator) provides obtaining a set of quasi-random sequences of nucleotides 
depending on user-defined input parameters (type of nucleic acid, sequence length, GC and dinucleotide composition). 
GATCGGenerator excludes the presence of any nucleotide repeats and homopolymer sites longer than three elements. The 
generated sequences can be used as service or masking sequences (e.g., in DNA steganography) and are suitable for any non- 
biological enzymatic manipulations. It is possible to generate numerous artificial nucleotide sequences and use them to create a 
universal oligotheca suitable for multiple encoding of non-biological data and their long-term storage. 

The data presented in Table 3 summarizes the results of the program. For a certain type of nucleic acid (in this case, 
DNA), the following data is shown: the content of dinucleotides (NN %), the number of generated sequences, their 
length (nucleotides — nt), and GC composition. 


Table 3 
Examples of short sequences differing in length, GC composition, and dinucleotide content (%) 
Input parameters Output parameters 
T Numb Length, | GC, | NN, Nucleotide sequence, 5’—3’ |Length,| GC, | NN, 
mber 
el {ont | % | nt | % | %* 


CTGG**TATATCGGAATCATATCGCGCAGTGT | 30 46.7 | 20.0 
AATCAGCTAGTAGGACGCAGTAGTGAATCA — {30 43.3 | 20.0 


a) a GAATGTAGTCCTAGGCACATACTACGTAGC = |30 46.7 | 20.0 
” AGTTGCACTGAAGTCTATGATCTGGCATGC | 30 46.7 | 20.0 

20 GACACACTACTATGGACGTGAGGCACTTAC  |30 50.0 | 20.0 
TCAGCTCAGCGCCAATCGAGCTTATAGTGC | 30 53.3 | 20.0 

51 GAGGCTATCGTCAAGCATAGACCGTGTGCT | |30 53.3 | 20.0 


5 30 60 GACTCAGTAGCTGCTCCGGACATACAGCCT  /|30 56.7 | 20.0 
TCGCGCGTTAGACTTAGGTCTCATCGCAGC | 30 56.7 | 20.0 
ACGCTCACAGGAGTTCGCATCGAACGATGC  |30 56.7 | 20.0 
ACGACAGTGATATAGCACGACGTGCTCATA | 30 46.7 |0.0 
GACTACATCTGATAGTACACGTGCTGCACT |30 46.7 |0.0 


DNA 5 a 0 | TCTATCTCTGCTAGAGCGCTCGTCACTCTA /30 50.0 |0.0 
TCTGATCTACTATAGCGATACGTGAGAGTG |30 43.3 | 0.0 
ACACATATATCGACGCACGCGTCGTAGTAC | 30 50.0 | 0.0 
TGCATGACCATGCTTGCGGTAGACATTCA 50 52.0 120.0 
GACGCGCGAATAGTAGGACGA 
GCATACGAGTGGCATACATATTAGACTAT 50 42.0 120.0 
ACGGTAGTGCATATGGTGCAA 

4l-— CTGAGACTCCTCTCTGTGGAGCTCCTAGTA 
: x 60 my CCGTCACGCGTGCTCTGAAG a cau ip 


CTGTGTGAACATACGATGCATTCTCATCTC 
GGTATGGCTGAAGTGCACAT 
GCGCTGACGTCATGGTTCATACCAATGTA 
GCATGATGTGCGATAGGCACA 

* NN shows the fraction (%) of the dinucleotides contained in the nucleotide sequence. 

** Dinucleotides are highlighted in bold. 


50 46.0 |20.0 


The obtained 30-nucleotide sequences were tested using the Blast tool from NCBI. The absence of 100% homology 
with known DNA sequences of living organisms was identified. The maximum coincidence was observed for the 
generated sequences with a length of 25 nucleotides (similarity is about 80 %). This indicates the ability of the 
GATCGGenerator to generate non-biological nucleotide sequences with high efficiency. It can be assumed that the 
sequences generated in this way do not have an absolute coincidence with the nucleotide fragments of living organisms. 
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In this case, special DNA-oligonucleotides of artificial origin containing informative and service parts can be used 
as a convenient information carrier. Recently, the authors of this work have proposed the use of NYRN- 
oligonucleotides [14] consisting of: 

— internal part (YR)n encoding the encrypted information; 

— service (auxiliary) parts S1 and S2 flanking sequence (YR)n (Fig. 2). 


> (N)x (YR)n (N)m 3? 
service encoding service 


part S1 part part S2 


Fig. 2. Structure of NYRN-oligonucleotides: N — degenerate nucleotides; Y — pyrimidines (C or T); R — purines (A or G); 
k, n, m — indices corresponding to the length of the part 


The length of the sites (n, k and m) may vary, but the structure of the service parts should provide the 
successful course of amplification reactions (length more than 18 nt, 40-60 % GC composition, absence of 
homopolymer sites and repeats). GATCGGenerator allows including NN dinucleotides containing identical paired 
nucleotides (e.g., AA, GG, CC, TT or UU for RNA), which can increase the specificity of molecular hybridization 
of nucleic acids. 

Discussion and Conclusion. Thus, based on the results of the scientific investigation performed, a software 
solution (GATCGGenerator) has been proposed, which, in comparison to traditional approaches, assumes more 
stringent conditions for generating sequences. Due to this feature, the limitations of known DNA generators are 
removed, and quasi-random sequences of nucleotides are formed depending on the specified input parameters. The 
obtained 30-nucleotide sequences were studied. The test allowed us to establish the absence of 100 % homology 
with known DNA sequences of living organisms. The generated sequences with a length of 25 nucleotides 
coincided as much as possible (by about 80 %). 

Note also that in order to hide information in NYRN-oligonucleotides, it is required to mix them with masking 
DNA. Masking sequences should be similar to sequences of NYRN-oligonucleotides, so that when trying to read 
hidden information, it would be impossible to recognize them without key sequences. The addressee should know 
the key sequences — primers to the service sites of NYRN-oligonucleotides. The addressee can decipher the 
transmitted message by isolating informative nucleotide sequences using a polymerase chain reaction followed by 
sequencing. A set of NYRN- and masking oligonucleotides can be easily obtained using GATCGGenerator, 
synthesized, and then stored as an oligotheca. To do this, it is enough to determine the optimal NYRN- 
oligonucleotides with subsequent filling of the oligotheca. In the future, it is planned to conduct laboratory 
experiments in order to test the proposed method of storing non-biological information and checking the viability 


of oligotheca obtained using the generator. 
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Abstract 

Introduction. Perianal fistula rapidly develops an abscess, requiring surgical decompression. However, simple cases must 
be managed. However, for patients with renal insufficiency, MRI with contrast is contraindicated. It is proposed to use 
diffusion-weighted images that can diagnose anal fistulae, showing areas of high signal intensity (inflammatory tissues). 
The aim is to determine sensitivity of diffusion-weighted image combined with T2 turbo inversion recovery magnitude 
and as an alternative technique to contrast-enhanced MRI using clinical examination as a reference. 

Materials and Methods. Study included fifty patients with a clinical diagnosis of perianal fistula. MRI sequences were 
T2 turbo inversion recovery magnitude in oblique coronal and axial planes, diffusion-weighted image, and T1 weighted 
image turbo spin echo (fat suppression) pre- and post-administration of contrast agents in oblique axial planes. Three 
radiologists evaluated the MR imaging data using a questionnaire of parameters that necessitated a binary response, “yes” 
or “no” answer. 

Results. Diffusion-weighted image combined with axial T2 turbo inversion recovery magnitude sequence had 96.7 %. 
All raters agreed that it is sensitive enough to correctly identify perianal fistula with a moderate Kappa agreement 
(k = 0.586) and p-value<0.001. The mean value of rater's responses was 76.7 % represents sensitivity of diffusion- 
weighted images + T2 turbo inversion recovery magnitude as an alternative technique to Tl-enhanced contrast with 
moderate (k = 0.553) agreement between raters and P-value<0.001. 

Discussion and Conclusion. Diffusion-weighted images and T2 turbo inversion recovery magnitude sequences exhibit 
comparable efficacy to T1-enhanced contrast sequences in detecting perianal fistula. This may be an option for patients 
with renal impairment who cannot receive an MRI contrast. 


Keywords: magnetic resonance imaging; diffusion-weighted image; T1-enhanced contrast 
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AHHOTalna 

Beeoenue. Ulapapextanbublit cBuu ObICTpO MepexogMT B aOcl{ecc, Tpebyromuli xupyprM4yeckoli Jekomupeccnu. Tem He 
MeHee, IIpoctTble cilyyau MosexaT JeueHHio. OgHaKO, WIA MalWweHTOB C MoYedHOM HeyoctTaTouHOcTEIO MPT c 
KOHTpacTHpoBaHvemM mpoTuBonoKa3zaHa. IIpeymaraetca mpuMeHeHue uddy3HOHHO-B3BeleHHOTO H300paxKeHHA, 
KOTOPoe CHOCOOHO HarHOCTHpoOBaThb aHasIbHble CBUUIM, WOKa3bIBad Ha CHUMKaX OOJIACTH C BbICOKOM HHTCHCHBHOCTBIO 
curHana (oy4ar BocnaneHua). Llemb paOoTbl 3akmIouaeTcA B TOM, YTOOBI OMpeyeUTb CTeIeHb UyBCTBUTeIbHOCTH 
nw300paxeHua Jupdpy3HOHHOH CieKTpaIbHOM TOMOrpaduH B COUeTAHHH C BEJIMYHHOM BOCCTaHOBJIeHHA TypOo-HHBepcHn 
T2 MW BO3MO2%KHOCT IIpHMeHeHHA aHHOrO MeTOJa B KayeCcTBe ajibTepHaTHBHOrO MeToya MPT c KoHTpacTHbIM 
yCHJICHHeM, HCMOIb3yA KIMHMYeCKoe OOcIeqOBaHHe B KavecTBe CTaHapTa cpaBHeHHA. 

Mamepuanoi u memooot. B xoge padotst Obi oOcreqoBaHbr 50 maljweHTOB C KMIMHMYeCKHM HarHo30M 
«napapeKTasIbHBI cBnup>. UocmeqopatembHoctu MPT mpeyctapianm coOol tTypOo-uHBepcuio T2 c BemM4HHOH 
BOCCTAHOBJICHHA B KOCbIX KOPOHapHBbIX HM aKCHAJIbHBIX IJIOCKOCTAX, H300paxeHue Zupdy3HoHHoi MPT u B3BelieHHoe 
w300paxeHHe ObICTpOro CIHH-9xa (0KHPOMOMABIIEHHe) NO H MOCIe BBEACHHA KOHTPaCTHBIX BEL[eCTB B KOCbIe aKCHaJIbHBIe 
TI0cKocTu. Tp paqMonora oueHuiu WaHHbre MPT, ucnomb3ya aHKeTy MapaMeTpoB, KOTOpad TpeOoBala OuHapHoro 
OTBETA, OTBETAa «a» HIM «HET». 

Pezyiomamoet ucciedosanua. Juddy3v0HHO-B3BelIeHHOe M300paxKeHHe B COUeTaHHH C MOCIIeOBaTeIbHOCTbIO 
BOCCTAHOBJICHHA aKCHaIbHOM TypOo-uHBepcun T2 moKazao0 96,7 %. Bce 9KcIepTbI COrsacHJIMCb C TeM, 4TO JaHHbIit 
MeTOJ, NOKa3aBLUMi yMepeHHoe cormacne Kanna (k = 0,586) u 3Hayenue p<0,001, qocraTouHo 4yBCTBUTeIeH WIA TOTO, 
YTOObI BEPHO HACHTHPUUMpOBATS MapapeKTaIbHbI cBuu. CpeqHee 3HavYeHHe IKCHePTHOM OLeHKH CocTaBHIO 76,7 %. 
OTO MOATBeEpxKLaeT, YTO YYBCTBUTeEJIbHOCTS ipesIO%KeHHOTO MeTOAa — Audpdy3HOHHO-B3BeLIeHHbIe U300paxKeHHA + 
BeJIMYHHa BOCCTaHOBJICHHA TypOo-HHBepcuu T2 — aABJIAeTCA aIbTepHATHBOM MeTOAY KOHTpacTHoro ycueHua T1 c 
yMepeHHBIM (k = 0,553) cormacveM Mexyy 9KcnepTamMu Hu 3HadeHHem P<0,001. 

O@6cystcdenue u 3akntouenue. JIuddy3MOHHO-B3BelICHHbIe M300paxKeHHA UM TOCeqOBaTebHOCTH BeJIMYMH 
BOCCTaHOBJICHHA IIpH TypOo-nHBepcun T2 eMOHCTpHpyIoT cpaBHUMy1o 9¢pdeKTHBHOCTh C MOCeAOBaTeIbHOCTAMU 
KOHTpacTHoro ycunenua Tl ya OOHapy2xXeHHA MapapektambHoro cBuuja. IpeqoxeHHbIi MeTOX MOXKeT CiLy2KHTb 
BapHaHTOM JIA WallMeHTOB C MOUeC4UHOH HEAOCTATOUHOCTBIO, KOTOPbIM MpoTHBONOKa3aHa MPT c KoHTpacTupoBaHHeM. 


KuroueBbie CI0Ba: MarHHTHO-pe30HaHcHad ToMorpadua, quddy3HOHHO-B3BellIeHHOe H300paxKeHHe, KOHTpacTHoe 
ycunenne T1 


Baarojapnocrn: aBTOpbI OnaroqapatT peakUWHt0 UH peleH3ecHTOB 3a BHHMAaTeCJIBHOe OTHOMCHHE K CTaTbe HU 3a YKa3aHHbIe 
3aMe4aHhs, MOSBOJIMBINIHE YIYAWINTb KaYeCTBO CTaTbH. 


Aaa wutTupopannusa. Noor Fadhil Baqir, Rasha Sabeeh Ahmed, Khalee Ibraheem Mohsen. Sensitivity of Diffusion- 
Weighted Image Combined withT2 Turbo Inversion Recovery Magnitude Sequence and as an Alternative to Contrast- 
Enhanced MRI in the Detection of Perianal Fistula. Advanced Engineering Research (Rostov-on-Don). 2023;23(3): 
307-316. https://doi.org/10.23947/2687-1653-2023-23-3-307-316 


Introduction. Perianal fistula is a rare condition that is difficult to treat and is associated with patient morbidity. It is 
the connection between the mucosal lining of the anal canal and the perianal epidermis [1]. The rate of occurrence is 
approximately two times higher in males than in females. The male-female ratio is 2:1. Perianal discharge that comes and 
goes, itchiness, discomfort, fever, and localized pain are the most typical symptoms [2]. Perianal fistula is associated with 
the rapid development of an abscess, necessitating timely surgical intervention for decompression. Consequently, the 
management of uncomplicated cases is crucial [3]. The treatment of anal fistula involves the excision of the primary 
fistula opening, any associated tracts, and any supplementary openings while ensuring the preservation of continence [4]. 
This requires accurate identification of the fistula's internal opening and any secondary abscesses or extensions. Imaging 
techniques, including magnetic resonance imaging (MRI), are increasingly important in perianal fistula diagnosis [5, 6] 
because they provide the most precise and reliable data regarding the anatomical characteristics of the fistula tract and its 
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association with the anal sphincter muscles. Providing these specific details by MRI has enhanced surgical success rates 
and reduced the likelihood of recurrence in perianal fistulas [7, 8]. T2-weighted with fat suppression (Turbo spin echo or 
turbo inversion recovery magnitude) MRI is very good at showing a fistula track. Still, some reports have pointed out 
how hard it is to distinguish between real abscesses that need to be drained by surgery and inflamed masses that cannot 
be drained and need anti-inflammatory medicine. Gadolinium chelate-enhanced T1-weighted MR scans are often used as 
a tool to help solve problems [9, 10]. Post-contrast MR imaging is the optimal sequence for detecting and characterizing 
perianal fistula. It is considered the gold standard by various institutions [5, 11]. 

Nonetheless, using gadolinium contrast agents increases expenses and duration and poses a risk of nephrogenic systemic 
fibrosis in patients already experiencing renal insufficiency. In addition, several studies have found evidence of gadolinium 
accumulation in the brain's deep nuclei after repeated contrast agent injections [1, 12, 13]. In recent years, Diffusion- 
Weighted Imaging (DWI) has provided high differentiation between pathological lesions, such as cancer or inflammatory 
processes, and the adjacent normal tissue [14]. DWI is a technique that captures variations in water mobility resulting from 
interactions with macromolecules, cell membranes, and changes in the tissue environment. Therefore, inflammatory tissues 
may also demonstrate limited diffusivity, resulting in the appearance of high signal areas on DWI, which takes advantage of 
this property. This sequence has the potential to be helpful in the diagnosis of anal fistulae [15, 16]. 

Moreover, Diffusion-Weighted Imaging can be integrated into customary Magnetic Resonance Imaging (MRI) scans 
of the perianal fistula, as it removes the need for a contrast agent, has shorter sequence duration, and does not necessitate 
additional equipment. The diffusion-weighted imaging technique has a naturally low spatial resolution, one of the 
method's drawbacks. Therefore, the T2 turbo inversion recovery magnitude (TIRM) sequence was added to DWI to assess 
the perianal fistula to overcome these drawbacks [16, 17]. In the T2 TIRM imaging technique, the fistulous tracts are 
commonly visualized as highly intense linear structures [18]. Additionally, it should be noted that pus, fluid, and 
granulation tissue manifest as areas of high signal intensity against a background of low signal intensity fat [6, 19]. 

Boruah et al. [20] evaluated the efficacy of diffusion-weighted magnetic resonance imaging (MRI) as a diagnostic 
tool for perianal fistulae. The study consisted of 47 participants diagnosed with perianal fistula. The MRI results were 
correlated with the clinical evaluation or the surgical findings. The performance of DWI alone was inferior to that of T2W 
imaging. Nevertheless, the most notable outcome was the assessment of combined DWI-T2W images, despite the lack of 
statistical significance compared to DWI or T2W images individually. 

Gu et al. [21] explored an alternative technique comparable to Gadolinium-contrast enhancement to assess any 
morphological changes in the preoperative evaluation of anal fistulas. The study comprised 46 individuals who had 
received a diagnosis of anal fistula. The gold standard for comparison will involve outcomes derived from surgical 
pathology. The utilization of FS T2-weighted and diffusion-weighted imaging exhibited similar effectiveness when 
compared to FS T1-weighted contrast-enhanced imaging for evaluating changes in the morphology of anal fistulas. 

This study aims to determine the sensitivity of diffusion-weighted images combined with T2 turbo inversion recovery 
magnitude (TIRM) and as an alternative technique to contrast-enhanced MRI using the clinical examination as a gold 
standard reference. 

Patients and Methods. The study was carried out between December 2022 and May 2023 at the MRI department of 
the oncology teaching hospital located in Medical City, Baghdad, Iraq. The research was carried out as a cross-section 
and utilized a non-randomized design. After a thorough clinical evaluation of each patient, the surgical outpatient clinics 
sent fifty patients who were thought to have perianal fistulas to the MRI department. Each participant gave their informed 
written agreement, and the institution's ethical committee gave its stamp of approval before the study was carried out. 

The inclusion criteria for this study include adult patients, regardless of gender, who have been diagnosed with a 
perianal fistula during a clinical examination and have been referred from outpatient surgical clinics. 

The exclusion criteria for this study include individuals who fall into specific categories, namely children, pregnant 
women, and patients who are contraindicated for MRI examinations. This includes patients with cardiac pacemakers or 
metallic shells and those who experience claustrophobia. 

The MRI examinations were performed using 1.5 T scanners (Magnetom Aera, Siemens Healthcare, Germany) 
equipped with a phased-array surface coil comprising 18 channels. There was no requirement for patient preparation 
before the examination. During the MRI, patients take a supine position on the examination table, with a phased array 
surface coil positioned over their pelvic region. It is crucial to ensure correct orientation when using magnetic resonance 
imaging (MRI) for imaging the anal canal. The planning of the images is performed with the assistance of a sagittal T2- 
weighted localizer that passes through the middle of the patient. The oblique transverse and coronal sections are 
orthogonally and parallel to the anal sphincter. 

The following MRI sequence parameters were utilized for this investigation: T2 TIRM in the coronal oblique plane 
(TR/TE: 5000/50, slice thickness: 3 mm, flip angle: 140°, interstice gap: 20 %, number of slices: 34, FOV: 390mm, matrix 
size: 384x288, signal averages: 1, inversion time (TI): 160 s, acquisition time: 2:27 (min: sec)). T2 TIRM performed in 
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the oblique axial plane (TR/TE: 5040/50 ms, slice thickness: 3 mm, interstice gap: 20 %, flip angle: 140°, number of 
slices: 36, FOV: 380 mm, matrix size: 384x234, signal averages: 1; inversion time (TI): 160 s, acquisition time: 2:02 
(min: sec). Single shot echo planar imaging (EPI) imaging (DWI) in an oblique axial plane (TR/TE: 2600/70, FOV: 
380mm, distance factor: 20 %, slice thickness: 3 mm number of slices: 30, matrix size: 128x128, b-value 50-400— 
800 with signal averages of | for 50, 3 for 400, and 800, acquisition time: 2:22 (min: sec)). And T1 FS TSE in an oblique 
axial plane before and after administration of contrast agent (TR/TE: 619/20, slice thickness: 3mm, number of slices: 27, 
interstice gap: 20%, flip angle: 160°, FOV: 239 mm, signal averages: 1, matrix size: 192163, acquisition time: 
1:55 (min: sec)). 

Three radiologists, each with an advanced radio diagnosis degree and over five years of experience in MRI reporting, 
were carefully selected for employment at the oncology teaching hospital based on their superior rankings among other 
radiologists. Furthermore, these individuals consistently produce a substantial volume of over 40 reports per week. They 
were provided information regarding the patient's clinical suspicion of perianal fistula while unaware of the surgical and 
prior MRI reports. The radiologists reviewed the MR imaging data sets using a questionnaire of parameters requiring a 
binary response (yes or no) as follows: 

— Diffusion-weighted image (3mm slice thickness) combined with T2 TIRM; 

— Diffusion-weighted image replaces T1 enhanced contrast with fat suppression. 

The data were inputted, categorized, and examined using the SPSS (statistical package for social sciences) software 
program version 26. The data were presented using frequencies and percentages. A significance level of 0.05 or lower 
was used to determine statistical significance. The study employed the Fleiss kappa test to evaluate the reliability and 
inter-rater agreement among multiple raters. 

Results. The findings of this study were derived from the analysis of a sample of fifty individuals who presented with 
a perianal fistula as determined through clinical assessment. Most individuals diagnosed with anal fistula were male, 
comprising 82 % of the cases, while females constituted the remaining 18 % as shown in Figure 1. The mean age of the 
entire patient population was 39.80 years, and the standard deviation was 12.46 years. 

According to the prepared questionnaire, the three rater's responses about the best visualization of perianal fistula by 
diffusion-weighted image combined with T2 TIRM and as an alternative method to T1 enhanced contrast. The frequency 
distribution of the three raters is demonstrated in Table 1. The three raters' responses according to the constructed 
questionnaire; regarding the best visualization of the perianal fistula compared to the positive clinical examination 
(standard reference) is represented in Fig. 2. 


=Female = Male 


Fig. 1. Frequency distribution of the gender for the 50 participant patients 
Table 1 
Frequency and percentage of DWI+T2 TIRM and DW1I as an alternative to T1 enhanced contrast images for three raters 


ides ahaa DWI image as an alternative to Tl enhanced contrast freq. (%) 
freq. (%) 
Rater 1 47(94 %) 37(74 %) 
Rater 2 49(98 %) 36(72 %) 
Rater 3 49(98 %) 42(84 %) 
Total 50 (100 %) 
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Fig. 2. Bar chart frequency distribution of the raters' opinions about using diffusion-weighted image combined with T2 TIRM and as 
an alternative to Tl enhanced contrast in detecting perianal fistula 


The average of the three raters for detecting cases of perianal fistula by each reported MRI sequence is demonstrated 
in Fig. 3. It was found that DWI+ axial T2 TIRM had a higher percentage (96.7 %), which corresponded to the sensitivity 
of this sequence in the detection of true cases of perianal fistula by the three raters. The average of raters' answers for 
replacing the T1 FS TSE enhanced contrast with DWI in detecting anal fistula was 76.7 % (sensitivity of DWI +T2 TIRM 
as an alternative technique to T1 enhanced contrast). The inter-rater kappa agreement was performed to test the reliability 
of the questionnaire used to evaluate MRI sequences obtained for the visualization of 50 patients clinically diagnosed 
with perianal fistula, as shown in Table 2. A significant (P-value<0.001) moderate (k=0.586) agreement was found for 
the question related to the use of optimized DWI+T2 TIRM. Also, significant (P-value<0.001) moderate (k=0.553) 
agreement between raters when asked to use DWI versus T1 enhanced contrast. 


m Average raters’ sensitivity 
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Fig. 3. Bar chart of overall average sensitivity of the raters’ opinions (%) 
about using different MRI sequences in detecting anal fistula 


Table 2 


Inter-rater reliability and overall agreement between radiologists in the diagnosis of anal fistula among 50 patients 


Sequence used Kappa 95% CI P-value 
DWI versus T1 enhanced contrast 0.553 0.548-0.558 <0.001 
DWI+ T2 TIRM 0.586 0.581-0.591 <0.001 


DWI diffusion-weighted image, TIRM turbo inversion recovery magnitude, CI confident interval. 


The presence of a perianal fistula can be seen in a variety of MRI sequences, as shown in Fig. 4 and Fig. 5, and Fig. 6 
demonstrates the abscess collection of patients with perianal fistula that shows the well-defined restricted diffusion of the 
abscess and hyper-signal intensity on T2 TIRM as an alternative technique to Tl enhanced contrast. 
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Discussion. The primary purpose of magnetic resonance imaging (MRI) in the context of perianal fistulas is to 
establish the extent of the tract, identify any side branches, and determine the presence of deep abscesses, especially at 
the supra levator level [22]. The utilization of Gadolintum-enhanced T1l-weighted images and fat suppression confers 
distinct advantages in differentiating a fluid-filled tract from an area of inflammation. The wall of the tract displays 
increased enhancement, whereas the central region exhibits a hypo-intense appearance. Post-contrast gadolinium images 
exhibit high efficacy in visually representing abscesses [23]. However, the gadolinium contrast agent does have certain 
limitations. These include an extended duration of magnetic resonance examination, additional costs, the potential 
deposition of gadolinium contrast in the brain, and contraindications for patients with renal insufficiency (1). 


c) 


Fig. 4. A 41-year-old female patient presenting with a left-sided supra sphincteric fistula at 2 o'clock (green arrow): 
a — shows the T1 FSTSE enhanced contrast; b — shows the T2 TIRM; 
c — shows the diffusion-weighted image (3 mm slice thickness); d — shows the ADC map 


Diffusion-weighted imaging has been proposed as a noteworthy alternative due to its ability to identify inflammatory 
lesions as areas exhibiting high signal intensity. The fistulous tracks demonstrate increased intensity in diffusion-weighted 
images, whereas the background signal is significantly attenuated. Consequently, the identification and extent of the fistula 
can be easily determined [24]. One limitation of diffusion-weighted imaging is its inherent drawback of having a relatively 
low spatial resolution. Diffusion-weighted imaging was not evaluated as an independent variable in our study. In this study, 
we evaluated the supplementary benefits of the technique in fat-suppressed T2-weighted imaging, specifically T2 TIRM 
(16). In our study, it was observed that all three raters concurred on the suitability of the majority of images derived from the 
sequences of Diffusion-Weighted Imaging (DWI) combined with T2 Turbo Inversion Recovery Magnitude (TIRM) to 
delineate anal fistula. Both rater two and rater three reported a percentage of 98, while rater 1 reported a slightly lower 
percentage of 94 for the images. Hence, the DWI+T2 TIRM sequence exhibited the highest percentage (96.7 %), aligning 
with the sensitivity of DWI combined with T2 TIRM compared to the clinical examination as the gold standard. The inter- 
rater kappa agreement was moderate (k = 0.586), with a significant P-value of less than 0.001. It is agreed Fahmy et al. (5) 
that is conducted a study that revealed a sensitivity of 92.12 % for the detection of perianal fistula when using a combination 
of diffusion-weighted imaging (DWI) and T2-weighted imaging (T2WI). 
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Fig. 5. A 51-year-old male presenting with a right-sided supra sphincteric branching fistula with extensive inflammatory changes and 
multiple extra sphincteric abscess collections (white arrow): a — optimized DWI; b — ADC map shows restrictd diffusion; 
c — T2 TIRM shows hyper signal intensity; d — T1FSTSE post-contrast shows peripheral enhancement. All the images are in the 
oblique axial plane 


Fig. 6. A 36-year-old male presenting with a left gluteal abscess collection with surrounding significant soft tissue edema (black 
arrow), a — diffusion weighted image; b — ADC map shows restricted diffusion; c— T1FSTSE post-contrast shows peripheral 
enhancement; d — T2 TIRM shows hyper signal intensity. All the images are in the oblique axial plane 


However, they regarded Tl-weighted imaging (T1WI) post-contrast as the reference standard for comparison. 
Additionally, Mohsen and Osman [25] evaluated the combined utilization of diffusion-weighted imaging (DWI) and T2- 
weighted (T2W) sequences to obtain optimal outcomes. The perianal fistulas and abscesses were accurately diagnosed in 
97.8 per cent of cases using T2-weighted imaging (T2WI) and T1-weighted imaging post-contrast (T1 WI post-contrast) 
as the established benchmark. Regarding the sensitivity of DWI combined with T2TIRM as an alternative technique to 
T1 FS TSE post-contrast, our research revealed that rater 3 said that it is sensitive for 42 (84 %), followed by rater one 
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and rater two in the same order, who reported 74 % and 72 % of the sensitivity, respectively. Therefore, an average raters' 
sensitivity of 76.7 %; there was a significant (P-value 0.001 fair (k=0.370) agreement between the three raters in selecting 
DWI as an alternative to post-contrast sequence for the detection of perianal fistula. Therefore, it reflects the capability 
of DWI in clearly delineating perianal fistula, making it comparable to contrast enhancement MRI in the detection of 
perianal fistula and possibly making it superior to contrast enhancement MRI in the detection of perianal fistula by 
reducing the amount of time an MRI examination takes and by reducing the adverse effects of contrast agents, particularly 
in patients with renal impairment. Our findings agree with those of Hori et al. (16), who established that there was not a 
statistically significant difference (p-value less than 0.05) in the sensitivity between combined contrast-enhanced 
sequences with T2WI and DWI coupled with T2WI. Also, according to the findings of Fahmy and Dawoud (5), the 
utilization of a combination of diffusion-weighted imaging (DWI) and turbo inversion recovery magnitude (TIRM) is 
deemed equivalent to the utilization of post-contrast images in the detection of primary and complicated fistula. The 
application of Diffusion-Weighted Imaging (DWI) has demonstrated its utility in differentiating between abscesses and 
inflammatory reactions and assessing the extent of disease activity. Active granulation tissue will exhibit enhancement 
on post-contrast Tl-weighted imaging compared to the fluid in the fistula. Rim enhancement indicates an abscess, while 
diffuse enhancement is typical of granulation tissue (5). In comparison, diffusion-weighted imaging is a technique that 
depicts changes in the mobility of water molecules caused by their interactions with cell membranes, macromolecules, 
and environmental changes in the tissue. On diffusion-weighted imaging (DWI), inflammatory tissues may exhibit 
reduced diffusivity, resulting in hyper-intense signal regions [26]. The apparent diffusion coefficient (ADC) values 
observed in abscesses are estimated to be lower than those in healthy tissues. As reported in the scientific literature, this 
is due to pus, which inhibits the unrestricted diffusion of water molecules (9). As inflammatory tissues typically 
demonstrate enhanced signal intensity, diffusion-weighted imaging may hold promise as a diagnostic sequence for anal 
fistulas (6). 

Our study has the limitation of not making a quantitative comparison by calculating ADC values in the diffusion- 
weighted images to diagnose active and inactive perianal fistula compared with the T1 enhanced contrast because of the 
overloaded work in the institution where the study was conducted. Therefore, we recommend future research to evaluate 
the ADC value of diffusion-weighted images to diagnose active and inactive perianal fistula compared with the Tl 
enhanced contrast with a larger sample size. 

Conclusion. The sensitivity of diffusion-weighted images in conjunction with T2 turbo inversion recovery magnitude 
(TIRM) sequences demonstrates a similar level of effectiveness as Tl enhanced contrast sequences in identifying and 
assessing perianal fistula. This combination of imaging techniques may serve as a viable alternative, particularly for 
individuals with renal impairment who are contraindicated to administer MRI contrast agents. 
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Abstract 

Introduction. The analysis of approaches to tracking the human body identified problems when capturing movements 
in a three-dimensional coordinate system. The prospects of motion capture systems based on computer vision are noted. 
In existing studies on markerless motion capture systems, positioning is considered only in two-dimensional space. 
Therefore, the research objective is to increase the accuracy of determining the coordinates of the human body in three- 
dimensional coordinates through developing a motion capture method based on computer vision and triangulation 
algorithms. 

Materials and Methods. A method of motion capture was presented, including calibration of several cameras and 
formalization of procedures for detecting a person in a frame using a convolutional neural network. Based on the 
skeletal points obtained from the neural network, a three-dimensional reconstruction of the human body model was 
carried out using various triangulation algorithms. 

Results. Experimental studies have been carried out comparing four triangulation algorithms: direct linear transfer, 
linear least squares method, L2 triangulation, and polynomial methods. The optimal triangulation algorithm 
(polynomial) was determined, providing an error of no more than 2.5 pixels or 1.67 centimeters. 

Discussion and Conclusion. The shortcomings of existing motion capture systems were revealed. The proposed 
method was aimed at improving the accuracy of motion capture in three-dimensional coordinates using computer 
vision. The results obtained were integrated into the human body positioning software in three-dimensional coordinates 


for use in virtual simulators, motion capture systems and remote monitoring. 
Keywords: motion capture, virtual reality, triangulation, computer vision, machine learning 


Acknowledgements: the authors would like to thank the Editorial board of the journal and the reviewer for their 


professional analysis and recommendations for correcting the text of the article. 


Funding information. The research was done on grant of the Russian Science Foundation No. 22-71- 
10057, https://rscf.ru/project/22-71-10057/ 


For citation. Obukhov AD, Dedov DL, Surkova EO, Korobova IL. 3D Human Motion Capture Method Based on 
Computer Vision. Advanced Engineering Research (Rostov-on-Don). 2023;23(3):317-328. 
https://doi.org/10.23947/2687-1653-2023-23-3-317-328 


© Obukhov AD, Dedov DL, Surkova EO, Korobova IL, 2023 


Information Technology, Computer Science and Management 


317 


http://vestnik-donstu.ru 


318 


Advanced Engineering Research (Rostov-on-Don). 2023;23(3):317—328. eISSN 2687-1653 


Hayunaa cmamoa 


Mertoy TpexMepHoOro 3axBaTa JBYKeHHH 4esIOBeKa HA OCHOBe KOMITIBIOTeCPHOrO 3peHnaA 


A.J. O6yxon D4, JL.JI. exon), E.O. Cypxosa\, H.JI. Kopodonal) 
TamOoBckHii rocyqapcTBeHHbI TeXHHYecKHM yHuBepcnuter, r. TaMOos, Poccniicxkaa Deyepalua 


DX] obuhov.art@gmail.com 


AnnoTayna 

Beedenue. poseyenupiii avasM3 CyWeCTBYIOWIMX MOAXOOB K OTCI@KHBaHHIO Tella YeOBeKa BbIABHI HaMune 
TipoOseM Ip 3axXBaTe JBMKeHH B TpexXMepHol cucTeMe KOOopsHHaT. OTMeyeHa MepCcieKTHBHOCTb CHCTeM 3aXxBaTa 
J{BWOKCHHM Ha OCHOBe KOMIMIbIOTepHOro 3peHHa. B cywecTByIOWJMx HCCIeOBaHHAX 10 OesMapKepHbIM CHCTeMaM 
3axBaTa JIBIDKeHHM paccMaTpHBaeTCA MO3HUMOHMpOBaHHe TOJbKO B JBYMepHOM IpoctpauctTBe. Ilosromy wWeJIbro 
HCCIeAOBAHHA ABJIAJIOCh MOBbILICHHe TOUHOCTH OMpeseweHuA KOOpAMHAaT ueNOBeyecKorO Tesla B TPeXMePHBIX 
KOOpMHaTax 3a C4eT pa3paOoTKH MeTOa 3axXBaTa J[BHXKCHHA Ha OCHOBeE KOMIBbIOTepHOrO 3peHHA HM asIrOPHTMOB 
TpHaHry AHH. 

Mamepuanoi u memooot. Ipenctapien MeToy 3aXBaTa JBWKCHHM, BKIOUAaIOWIM KasIMOpOBKy HeCKOJIbKUX KaMep U 
cbopMasM3alHIo Mpouesyp oOHnapyxXeHHA 4eOBeKa B KajIpe C UMCHOML30BaHHeM CBepTONHON HelipoHHon ceTu. Ha 
OCHOBE IIOJIYYCHHBIX OT HeMPpOHHOM CeTH CKEJICTHBIX TOUCK OCYLIECTBIIACTCA TPeXMepHad PCKOHCTPyKUMA MOJeIM TesIa 
YeIOBeKa C HCHOJb3OBAHHEM pa3JIM4HBIX aIrOpHTMOB TpHaHTry AHH. 

Pe3ynavmamot uccredoeanua. \[popeyeub! 9KCIepHMeHTAaIbHbIle HCCIeqOBaHHA TO CpaBHeHHIO YeTbIPex asIrOpHTMOB 
TPHaHTyJAWAM: MpAMoro JIMHeMHOrO MepeHoca, JMHeMHOrO MeTOa HaMMeHBbIUMX KBaypaTos, L2 Tpuanrynauuu vu 
TIOIMHOMHaIbHOrO) §=9MeTOZ0B. OnpeyeseH ONTHMaJIbHbI asITOpHTM TpHaHryIAWMH (MOJMHOMMAJIbHBIM), 
oOecreunBaloluii WorpeliHocTs He Oosee 2,5 nukceset uu 1,67 caHTuMeTpoB. 

O6cysicoenue U 3aKiIoveHUe. BlaBeHbl HeOCTATKM CyIeCTBYIOWIMX CHCTeM 3axBaTa JBWKeHHA. IIpeqnoxKeHHblit 
MeTO HalipaBsIeH Ha MOBbILICHHe TOYHOCTH 3aXBaTa J[BHXKeHHi B TPeCXMePHbIX KOOPAMHaTaxX C HCMOJIb30BaHHeM 
KOMIIbIoTepHoro 3peHHa. TlomydeHHble pe3yIbTaTbI MHTErPHPOBAHBI B IporpaMMHoe obecriereHHe NO3HIMOHMpOBaHHA 
Tela YeOBeKa B TPCXMEPHBIX KOOPAMHaTax JIA yaIeHHOrO MOHHTOPHHTa, UCHONb30BaHHA B BUPTYAJIbHBIX 


TpeHakepax WU CHCTeCMax 3axBaTa TBWWKeHMI. 


Korouesble cil0Ba: 3aXBaT TBIWKeHHH, BUpTyasIbHadA PeCayIBHOCTb, TPpHaHTyJIAWHA, KOMITbIOTepHOe 3PeCHHA, MaliHHHOoe 
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Introduction. Significant progress has currently been made in the domain of computer vision. Technologies have 
been developed to solve the problems of detecting objects, determining their state, geometric evaluation of the space 
depicted on the frame, and a lot more. As a result, computer vision has become widespread in various spheres of human 


activity, ranging from healthcare and education to entertainment. A rather promising direction is the use of computer 
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vision technologies for three-dimensional reconstruction and positioning of various objects, including people. There is 
fairly large number of systems for determining the absolute position of a person in space, which can be divided into the 
following categories. 

— Systems using inertial sensors and providing the determination of the amount of their movement, as well as the 
change of angles between them, which involves the use of gyroscopes and accelerometers [1]. A well-known 
representative of this category is the Noitom Mocap Perception Neuron [2], which includes up to 32 inertial sensors. 

— Laser position tracking systems based on the use of base stations installed on opposite sides of the room and 
emitting infrared rays, which provide accurate determination of the position and orientation of sensors in space. An 
example of such systems is a virtual reality kit from HTC [3], which have an error of up to 0.1 mm. 

— Systems using magnetic sensors [4] based on the use of a magnetic field to capture human movement, which 
assume the presence of wearable sensors on the user's body. This category includes Polhemus Liberty — a portable 
electromagnetic motion tracking system, considered one of the fastest (sampling rate — 240 Hz). 

— Marker-based optical systems determine the position of objects by markers using a set of cameras. An example is 
Vicon, which has a fairly low error: the average absolute errors of marker tracking are 0.15 mm in static tests, and 
0.2 mm (with corresponding angular errors of 0.3°) in dynamic tests [5]. 

— Marker-free optical systems based on the use of computer vision and machine learning. Examples of such 
technologies are Open space, MediaPipe, Movenext. With their help, human movements can be tracked with an 
accuracy of up to 30 mm [6]. 

After analyzing the listed categories of motion capture systems, it can be concluded that most of the solutions used 
to recognize human actions and movements involve various wearable devices, such as sensors or gloves. Most of these 
devices are bulky due to the large number of sensors and the need for a wired connection. Some systems have high 
accuracy, but they cannot be used due to the size or the presence of electromagnetic interference [7]. Inertial systems 
have a number of problems associated with the accumulation of errors, which limits their use only to relative 
positioning in space. 

Therefore, optical systems for recognizing and tracking user actions are well regarded. To get information about the 
actions and position of the user, frames obtained from the camera are used. Among optical systems, it is worth noting 
those that use markers (the user may be wearing special clothes or certain labels fixed on him), which makes it difficult 
to use them under real conditions. They are more applicable to specially prepared premises (e.g., film studios). 

Systems that do not use any markers allow users to interact more freely with the environment and are more suitable 
for use under real conditions. The significant disadvantages of systems in this line include relatively low accuracy, 
unreliability, and low performance. To a great extent, this may be due to the shortcomings of computer vision 
algorithms used to recognize a person in the frame, as well as the following reasons: the variability of a person’s 
appearance and lighting conditions, partial occlusions owing to the layering of objects in the scene, the complexity of 
the human skeletal structure. 

As a tule, the operation of marker-free motion capture systems is based on an algorithm for evaluating a person's 
posture. Approaches to solving the problem of assessing a person's posture can be divided into top-bottom and bottom- 
up. In top-bottom approaches, first there is a detection of people in the frame, then an assessment of the pose of each 
person found. Algorithms that relate to the bottom-up approach, at the first stage, search for body parts in the frame, 
then group them into poses. As a rule, convolutional neural networks are used for this task, such as YOLO (You Look 
Only Once) [8], SSD (Single Shot Detection) [9], R-CNN (Region CNN) [10], and others. They provide the recognition 
of numerous different objects, including a person or individual body parts with high accuracy. However, one of the 


disadvantages of the solutions listed above is their low performance and slow operation. To solve this problem, there 
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are special frameworks (MoveNet [11], MediaPipe [12], OpenPose [13]) that also use neural networks optimized for 
real-time operation. 

It should be noted that the above algorithms, technologies and approaches of marker-free motion capture systems 
provide positioning in two-dimensional space, which makes it difficult both to determine the distance to objects and 
their sizes, and to track complex movements when, e.g., the user's hands are hidden by his body. Existing solutions in 
the field of stereo cameras can be effective, but they are not very accurate when the object is significantly removed from 
the camera, which happens when tracking the entire human body. In addition, they do not solve the problem of 
occlusions. Thus, the major line of research is the development of a method of motion capture using multiple cameras 
and computer vision technologies. When implementing multi-camera motion capture systems, the problem of 
combining objects from several images inevitably arises, i.e, the need to perform triangulation. Among the 
triangulation methods, linear and iterative linear algorithms can be distinguished. 

Linear triangulation is the most common approach to performing reconstruction of objects in three-dimensional 
space, including such methods as linear-proprietary method, linear least squares method, direct linear transformation, 
which differ in varying degrees of resistance to noise [14]. 

Iterative linear methods are a more robust version of linear triangulation. Conventional linear methods may be less 
accurate when solving problems of triangulation of a set of points, since in this case, the minimized error has no 
geometric meaning (it does not take into account the shape of the skeleton and the rules for connecting points). The key 
idea of iterative linear methods is to adaptively change the weights of linear equations in such a way that the weighted 
equations correspond to errors. Iterative linear methods include L2 and Loo triangulation [15]. 

Thus, within the framework of this study, the following task was set: to develop a method for capturing human 
movements that provides positioning the user's body in three-dimensional coordinates with minimal error and using 
computer vision technologies. The proposed method can be used as a replacement for existing motion capture systems, 
or as part of other algorithms, e.g., for the subsequent classification of a person's condition. This work was aimed at 
increasing the accuracy of determining the poses and coordinates of the human body in three-dimensional coordinates 
by developing motion capture methods based on computer vision. To achieve this goal, it was required to formalize the 
main stages of the process of capturing points of the human body from several cameras, integrate triangulation 
algorithms, choosing among them the optimal one from the point of view of accuracy, carry out the software 
implementation of the proposed method. 

Materials and Methods. Solving the problem of 3D positioning of a person in space includes the following main 
stages: 

— preliminary calibration of a set of cameras; 

— implementation of human detection procedures in the frame, and calculation of skeletal points; 

— calculation of 3D reconstruction of the human body model. 

Let us look at them in more detail. 

The calibration process involves the camera system taking several pictures of a calibration template, on which it is 
easy to identify key points with known relative positions in space. After that, internal and external parameters are 
calculated for each camera. Internal parameters are constant for a particular camera, external parameters depend on the 
location of the cameras relative to each other [16]. Therefore, this step must be performed before the first use of the 
camera system in a given location. 

To calculate the coordinates of a point in three-dimensional space, it is necessary to know the coordinates of its 
projection on the images and the projective matrices of the cameras [17]. Projective matrix P of some camera can be 


represented as a combination of matrices A (containing the internal parameters of the camera) and R (rotation), as 
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well as the displacement vector T , which describe the change of coordinates from the world coordinate system to the 


coordinate system relative to the camera: 


f. 0 c.f te Ns 4 
P=AR|T]=| 0 Fy || ty Me. Be be ls (1) 
0 0 ho Te Ty G 


where (x, y) — coordinates of the projection of a 3D point on the image in pixels; (c,,c,) — coordinates of the central 
point of the camera; (f,, f,,) — focal length in pixels. 


At the second stage, it is required to obtain directly the key (skeletal) points of the human body on each of their 
cameras. To extract skeletal body points from the frame, it is possible to use various machine learning technologies, 
e.g., MoveNet, MediaPipe, OpenPose, and others [18]. As part of this study, it is proposed to use a highly efficient and 
productive Pose module from the MediaPipe library. MediaPipe Pose uses machine learning to accurately track a 
person's body posture, determine 3D landmarks, and mask background segmentation on the entire body from RGB 
video frames. This approach makes it possible to track up to 33 points and provides real-time operation on most modern 
devices. 


Thus, as part of the second stage, a set of 33 points is formed for each i -th camera: 
{xy = (uy Vy) | 7 ML 25--033},7 € (1,2... KH, (2) 
where u,,— coordinate of j — th point on X axis in i-th image; v, — coordinate of j-th point on Y axis in i-th 


image; K — total number of cameras and images. 

At the third stage, the positions of key skeletal points in three-dimensional space are calculated. To obtain data on 
the position of human skeletal points in space, triangulation is performed — finding the coordinates of a 3D point by the 
coordinates of its projections. Triangulation is one of the most important challenges in computer vision, its solution is a 
crucial stage in 3D reconstruction, it affects the accuracy of the entire result [19]. 

Epipolar geometry is fundamental for the three-dimensional reconstruction of the object points based on the position 
values of the projections of the points in the images from all cameras. Its main idea is that 3D points in the scene are 
projected onto lines in the image plane of each camera — epipolar lines. These lines correspond to the intersection of 
the image plane and the plane passing through the camera centers and the 3D point. This idea provides a condition for 
finding pairs of corresponding points on two images: if it is known that point x on the plane of the first image 
corresponds to point x' on the plane of another image, then its projection should lie on the corresponding epipolar line. 
According to this condition, the following relation will be valid for all corresponding pairs of points x <> x': 

x'Fx=0, (3) 
where /’ — fundamental matrix having size 3x3 and rank equal to 2. 

For some point X , given in three-dimensional space, the following projection formula expressed in homogeneous 
coordinates is valid: 

Px, (4) 
where x, = w(u,,v,,1)7 — homogeneous coordinates of some point on the plane of the i-th image (obtained from the 


i-th camera during the second stage), including the position on image u, (on X axis) and v, (on Y axis); w — scale 


factor; P — projection matrix of i -th camera obtained at the first stage. 
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To simplify calculations, the projection matrix of the camera is often presented in the following form: 


P =| pt |(PeR*), (5) 


where p/’ — j-th row of matrix P. 


Therefore, equation (4) can be represented as follows: 


wu, = pi X, 
wy, = pet x, (6) 
we p37 X. 


Since w — scale factor, we obtain the following system of equations: 


u,pi"X ~ pit X =0, , 
u, p37 X — pr X =0. 4 


Since X is a homogeneous representation of coordinates in three-dimensional space, then, for their calculation, it is 


necessary to obtain x, and P for at least two cameras. To solve the system of equations (7), 4 algorithms were 


considered [14]: 
— direct linear transfer (DLT); 
— linear least squares method; 
— L2 triangulation; 
— optimal (polynomial) method. 
DLT refers to a linear triangulation algorithm, whose main advantage is the simplicity of its implementation. 


Specifically, in the OpenCV computer vision library there is a ready-made implementation of this algorithm in the 


triangulatePoints method. 


The linear least squares method also refers to linear ones and consists in the fact that the system of homogeneous 
equations (7) is reduced to a system consisting of inhomogeneous equations, for whose solution, the least squares 


method is used. 


L2 triangulation is an iterative method of three-dimensional reconstruction, whose solution is reduced to minimizing 
the reprojection error: 


Yd(x,,x,) > min, (8) 


where x, — coordinate of the projection of the estimated point in the image; x; — projection coordinate calculated 


from formula (4) for an already determined spatial point; d(e) — distance between two points. 


The algorithm of optimal (polynomial) triangulation refers to non-iterative approaches. To solve it, a sextic 
polynomial is required. The minimization criterion for performing three-dimensional reconstruction in this method can 


be defined as follows: 
Ed(x,,2,) > min, (9) 


where 2, — epipolar line corresponding to point x, . 


When using a two-camera system, to minimize error (9), the following sequence of actions must be performed: 


— parametrize the bundle of epipolar lines in the first image using parameter ¢. Thus, the epipolar line in the first 
image can be expressed as 1, (f) ; 
— using fundamental matrix /F’, calculate the corresponding epipolar line 1,(¢) in the second image; 


— express the distance function (9) as a function of f ; 


— perform a search for value ft, at which (9) tends to a minimum. 
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Using the methods of elementary calculus, it is possible to reduce the solution of the minimization problem to 
finding the roots of a sextic polynomial. The calculation of the assumed spatial point is performed using the direct linear 
transfer method (DLT) [17]. 

Summing up the third stage, we get that after successfully solving system (7) and obtaining the world coordinates of 
the key points of the target object (human body), the following set of points H is formed: 

H ={X,|vi(x,=PX,)} (10) 
where X, — world coordinates of the skeletal point of the human body obtained after solving the triangulation 
problem, expressed in centimeters. 

Thus, in this study, the optimization problem, when using two cameras, is reduced to finding triangulation method 


MT : {x,} > H., in which the reprojection error tends to a minimum: 


R=2 > min. (11) 


Research Results. Optimization problem (11) is solved through performing triangulation of 2D object points 
obtained from images of several cameras, in the framework of this study — from two cameras using various algorithms 
listed in the previous section. 

The listed triangulation methods were implemented using OpenCV and NumPy libraries. For comparison, the 
algorithms were integrated into software implementing the method of three-dimensional motion capture. An example of 


the method for reconstructing the entire human skeleton is shown in Figure 1. 


First camera 3D skeleton ff Second camera 
100 


Fig. 1. Example of the method, including recognition of a person on two cameras and construction of a 3D skeleton 


Then, these algorithms were compared by the value of the reprojection error function (11) for all points of the 
skeleton from two images. The comparison of the selected triangulation methods by the error rate, as well as by the time 
of obtaining a solution (computational complexity) for the entire set of skeleton points was carried out. Summary 


comparative diagrams are shown in Figure 2. 
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Fig. 2. Comparison of triangulation methods by metrics: a — by reprojection error; b — by calculation time 


A number of experimental tests were also carried out for the selected triangulation methods. Under testing, the 
calculated lengths of the user's limbs and the absolute deviation of the obtained values from the real ones were 
measured for each approach. The comparison is presented in Table 1. 


Table 1 
Comparison of the accuracy of determining the size of limbs in the process of triangulation 
Body segment DLT Least Squares L2 Polynomial Real value 
Forearm 25.2+ 1.6 30.8 + 0.2 26.6 + 0.5 24.3+0.4 26 
Shin 42.2+2.0 65.34 1.1 44.6+0.7 38.7 + 1.8 41 
Hip 45.7+£2.7 59.5+0.49 48.7 + 1.3 44.1+0.6 45 
Average deviation 2.43 14.58 2.26 1.67 0 
Presented are the average values (in centimeters) after a sample of 10 measurements + standard deviation in the sample 


The developed software includes the following modules: 


— for working with input devices (cameras); 

— to perform calibration and obtain basic camera parameters; 

— to synchronize multiple cameras; 

— for object recognition (user's body and arms); 

— to analyze the location of the found skeletal points; 

— to build real-time visualization. 

When implementing the software, the Python programming language, OpenCV and Matplotlib libraries were used. 
The operation of the system was carried out in several streams: one was responsible for receiving data from cameras, 
the second — for visualization, the third — for sending the received world coordinates of the human body to external 
systems or modules. Using a unified protocol with a data package in JSON format provides integrating the software into 
third-party systems (e.g., Unity game development environments, Unreal Engine, etc.) [20, 21]. 

Discussion and Conclusion. Let us analyze the results of comparing triangulation algorithms by selected metrics, 


shown in Figure 2 and in Table 1. 
During the comparison, it was found that the optimal algorithm for three-dimensional reconstruction was the 
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was no more than 3 %, taking into account the fact that MediaPipe Pose did not fix the upper point of the head and it 
was calculated approximately based on the position of the eyes. When measuring limbs, the error ranged from 0.9 cm to 
2.3 cm, the average was 1.67 (Table 1). Thus, real tests validate the correctness of the choice of the polynomial method. 

Next, we compared the results obtained with existing studies, e.g., described in [22]. The authors also used trained 
networks (OpenPose) to implement a marker-free human recognition system, a camera calibration procedure, and the 
extraction of skeletal points, but placed cameras next to each other to simulate stereo vision. This key difference made it 
possible to recognize human postures within the framework of this study, when some parts of the body overlapped 
others. In addition, using MediaPipe Pose provided tracking 33 skeletal points, not 18, as in the OpenPose-based 
method. The obtained error values generally corresponded to existing studies (the best result in [22] was 2 cm), which 
allowed us to conclude that the proposed approach can be used in practice. Other marker-free systems, e.g., based on 
Kinect [23], also showed comparable results in terms of measurement error (2-5 cm). Thus, the resulting solution 
generally corresponded to the accuracy of existing developments. 

A comparison of the calculation time of a set of points, shown in Figure 2 on the right, demonstrated that the DLT 
algorithm provided the highest performance. However, all algorithms showed acceptable results (to provide a speed of 
30 and even 60 frames per second). Therefore, this metric was not determinative. 

The developed software can be used in various subject areas primarily as a replacement for motion capture systems 
based on inertial sensors. The advantages of the proposed solution are low economic costs for implementation and 
accessibility (transition from highly specialized motion capture suits to common camera-based tools), the possibility of 
parallel capture of body models of several users [24]. 

The scientific novelty of the research consists in a comprehensive approach to formalizing the process of three- 
dimensional positioning of a person using computer vision technologies. It includes preliminary calibration of a set of 
several cameras, formalization of procedures for detecting a person in a frame using an arbitrary neural network to 
obtain skeletal points, as well as calculation of three-dimensional reconstruction of a human body model using various 
triangulation algorithms. The study presents all the necessary calculation formulas and detailed steps to achieve the 
goal — to increase the accuracy of determining the poses and coordinates of the human body in three-dimensional 
coordinates using computer vision technologies. The theoretical results obtained are quite universal and can be used for 
the practical implementation of motion capture systems based on various models of neural networks, and not just 


MediaPipe Pose. 
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Abstract 

Introduction. Environmental problems arising in shallow waters and caused by both natural and man-made factors 
annually do significant damage to aquatic systems and coastal territories. It is possible to identify these problems in a 
timely manner, as well as ways to eliminate them, using modern computing systems. But earlier studies have shown 
that the resources of computing systems using only a central processor are not enough to solve large scientific problems, 
in particular, to predict major environmental accidents, assess the damage caused by them, and determine the 
possibilities of their elimination. For these purposes, it is proposed to use models of the computing system and 
decomposition of the computational domain to develop an algorithm for parallel-pipeline calculations. The research 
objective was to create a model of a parallel-conveyor computational process for solving a system of grid equations by a 
modified alternating-triangular iterative method using the decomposition of a three-dimensional uniform computational 
grid that takes into account technical characteristics of the equipment used for calculations. 

Materials and Methods. Mathematical models of the computer system and computational grid were developed. The 
decomposition model of the computational domain was made taking into account the characteristics of a heterogeneous 
system. A parallel-pipeline method for solving a system of grid equations by a modified alternating-triangular iterative 
method was proposed. 

Results. A program was written in the CUDA C language that implemented a parallel-pipeline method for solving a 
system of grid equations by a modified alternating-triangular iterative method. The experiments performed showed that 
with an increase in the number of threads, the computation time decreased, and when decomposing the computational 
grid, it was rational to split into fragments along coordinate z by a value not exceeding 10. The results of the 
experiments proved the efficiency of the developed parallel-pipeline method. 

Discussion and Conclusion. As a result of the research, a model of a parallel-pipeline computing process was 
developed using the example of one of the most time-consuming stages of solving a system of grid equations by a 
modified alternating-triangular iterative method. Its construction was based on decomposition models of a three- 
dimensional uniform computational grid, which took into account the technical characteristics of the equipment used in 
the calculations. This program can provide you for the acceleration of the calculation process and even loading of 
program flows in time. The conducted numerical experiments validated the mathematical model of decomposition of 
the computational domain. 
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AHHOTauHA 

Beedenue. Ikonornyeckue TpoOMeMbl, BOZSHUKAIOWIMe Ha MEJIKOBOJHBIX BOJOEMAX H BbI3bIBACMBIC KaK IIPHPOAHbIMH, 
TaK MW TeXHOreHHbIMM (bakTOpaMH, exKerOqHO HaHOCAT CyIeCTBeHHbIM yulepO akBacHcTeMaM uM UpHOpexkHbIM 
TeppuTOpHAM. CBoeBpeMeHHO OMpeeUTb ITH MpOOIeMbl, a TaKKe MYTH UX yCTpaHeHHA BO3MO2%KHO C MCHOIb30BaHHeM 
COBPeMCHHBIX BbIYMCJIMTeIbHBIX cHcTem. Ho mpoBeyzéHHble paHee UccleqOBaHHA WOKa3aIM, 4TO pecypcoB 
BBIYHCIIMTCIBHBIX CHCTeM, MCIOJb3YIOWIMX TOJIbKO WeHTpaIbHbIN Mpoweccop, HeAOcTaTOUHO WIA pewieHHA OOSIIHX 
Hay4HBIX 3a]{a4, B YACTHOCTH, NO NPOrHOSHPOBaHHIO KPYIHBIX IKOJOTHYCCKHX MPOMCLIeCTBHU, OWeHKe HaHeCceHHOrO 
MMU yijepOa HM OpeseseHuIO BOSMOX%KHOCTeH UX ycTpaHeHHaA. Ja 9TUX Wee NpelaraeTca MCHOML30BaTb MOEN 
BBIYHCIMTeIbHOM CHCTeEMbI HM JeKOMNO3HUMH pacuéTHOH oOsacTH WIA pa3spaOoTKM asIrOpHTMa MapaJiseIbHO- 
KOHBeHepHBIX BbIHCHeHuH. Lenbro WaHHOM paOoTEI ABJIAeTCA CO3qaHMe MOJeIM MapaliesIbHO-KOHBelepHoro 
BBIYHCIMTeIbHOTO Uporecca JI PeWICHHA CHCTeMbI CeTOUHBIX ypaBHeHHi MOAMPMUMpOBAaHHbIM TOMepeMeHHo- 
TPCYyIOJIbHbIM HTepallMOHHbIM MeTOJOM C HCHOb30OBAHHeEM J[CKOMMO3HIMH TpéxMepHo paBHOMepHoH pac4éTHOH 
CeTKH, YUHTHIBaIOLNeH TEXHHYECKHe XapakKTePHCTHKH UCMOIb3yeMOro JIA pacueTOB OOopyAOBaHHA. 

Mamepuanoit u memoovi. PaspaOotaHbl MaTeMaTHYecKHe MOJ{eIH BBIYMCIMTeIbHOM CHCTeMbI HM pacuéTHOM CeTKH. 
Mogetb eKOMMO3HUHM pacyéTHOM oOsacTH BbINOMHeHa C YyYETOM XapaKTepHCTHK TeTepOreHHOM CHCTeMBI. 
IIpequoxKeH MapasiesbHO-KOHBeHepHbIM MeTO PpelleHHA CHCTeMbI CeTOUHbIX ypaBHeHHi MOMHUMpOBaHHbIM 
HOMePeCMeHHO-TpeyFOsIbHbIM HTepallHOHHbIM MeTOOM. 

Pezyibmamei uccnedoesanua. Ha a3n1ke CUDA C nanmcana nporpamMa, peasmM3yrollad MapasviesIbHO-KOHBeHepHbIit 
MeTO PeLIeHHA CHCTeMbI CeETOUHBIX YpaBHeHHit MOAMPUUMpOBaHHBIM MOMePCMeCHHO-TpeyTOJbHbIM HTepal|MOHHbIM 
MeToyoM. IIpopewéHHbIe SKCIIepHMeHTbI NMoOKa3aIu, 4TO C YBeJIMYCHHeEM UHCa MOTOKOB BpeMA BbIMHCIeHHM 
yMeHbIUaeTCA HM TPH JeCKOMMO3HIMH pacuéTHOH CeTKH pal[MOHaJIbHbIM ABJIAeTCA pa30HeHHe Ha (pparMeHTBI 10 
KoopyuHate Z Ha BeIM4MHYy, He MpeBbiiuarouryio 10. PesynbtTaTbl 9KCIePHMeHTOB MOATBepAMIH 3PPeKTHBHOCTE 
pa3paOoTaHHOro MapasiesIbHO-KOHBeMepHoro MeTOsAa. 

O@6cyscdenue u 3aknio4uenua. Ilo uToraM IMpoBeyeHHbIX MCCIeOBaHHit pa3spadoTaHa MOJeyIb MapasiebHO- 
KOHBeMepHOrO BBIYHCIIMTeIbHOrO polecca Ha MpHMepe OAHOTO W3 CaMbIX TpyOEMKUX 3TANOB PeWICHHA CHCTeMBI 
CeTOYHBIX ypaBHeHHM MOAMPUUMpOBaHHBIM TOMepeMeHHO-TpeyroJIbHbIM HTepallMOHHbIM MeTOOM. Eé moctTpoeHue 
OCHOBaHO Ha MOJeIAX J©KOMNO3HUMH TpéxMepHOH paBHOMepHOH pacuéTHOM CeTKH, YYNTIBaIOLWelH TeXHHyecKHe 
XapaKTe€PpHCTHKH HCHOIb3yeMOro B pacueTax oOOopynoBaHuaA. IIpHMeHeHve IporpaMMBI MO3BOJIMT YCKOPHTb Tporecc 
pacuéTa WM paBHOMepHO 10 BpeMeHH 3arpy3HTb MporpaMMuble MoTOKH. I[poBexeHHbIe 4YHCIICHHbIe IKCIIEPHMCHTHI 
MOATBEP AHI MaTeMaTHYeCKyIO MOJeIIb JEKOMMO3HIMH pacuéTHO oONacTH. 


KoroueBble CJ10Ba: TlapasWIeJIbHBIM asIrOpHTM, BbIUMCIMTeIbHbIM Tipowecc, CCTOUHbIC ypaBHeCHHA 


BsaarojapHoctn: aBTOpbIl BbIpaxKatoT OlaroyjapHocTb peqakWHoHHOH KOJIICrHH WKYypHasiIa HU PeileH3eHTy 3a 
IIpod@eccvoHasIbHBIii aHasIn3 HW peKOMeCHaunu JIA KOPppeCKTHPOBKH CTaTbH. 


@uHancupospanne. Padota BEInouHeHa IpH MoszWepxKe Pocculickoro Hay4dHoro (ousa (mpoext Ne 21—71—20050). 


Aaa uuTuposanna. JIutaunos B.H., Pyneuxo H.b., [pauesa H.H. Paspadbotka MoyesIM MapasiesIbHO-KOHBeMepHOro 
BBIYHCIMTeIbHOTO Ipolecca Aid PeUWIeHHA CHCTeMbI CCETOUHBIX ypaBHeHui. Advanced Engineering Research (Rostov- 
on-Don). 2023;23(3):329-339. https://doi.org/10.23947/2687-1653-2023-23-3-329-339 


Introduction. Recently, a number of serious environmental problems have been observed in the Rostov region. 
These include, in particular, the eutrophication of waters of the Sea of Azov and the Tsimlyansk reservoir, which causes 
the growth of harmful and toxic species of phytoplankton populations [1]. Engineering works in the waters of rivers and 
seas cause pollution of adjacent territories, changes in the population structure of biota, and deterioration of 
reproduction conditions of valuable and commercial fish. Climate change in the south of Russia has led to an increase in 
the number of cases of flooding of some territories in the area of the Taganrog Bay and the floodplain of the Don River 
caused by up and down surges. In the last decade, during the summer period, almost complete drainage of the Don 
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riverbed was observed several times, which led to a complete stop of navigation. To predict the occurrence and 
development of such cases, to plan ways to address their consequences, to assess the damage caused by them, modern 
software systems built using high-precision mathematical models, numerical methods, algorithms and data structures 
are needed [2]. 

Mathematical models used in predicting natural and man-made disasters are based on systems of partial differential 
equations, such as Poisson, Navier-Stokes, diffusion-convection-reaction, and thermal conductivity equations. The 
numerical solution to such systems causes the need for operational storage of large amounts of data (in arrays of various 
structures) and the solution to systems of grid equations of high dimension exceeding 10°. The amount of RAM 
required to store arrays of data when numerically solving only one Poisson equation for a three-dimensional domain 
with a dimension of 103x103x103 by an alternately triangular iterative method is more than 64 GB. In the case of 
numerical solution to combined tasks, hundreds of gigabytes of RAM are required, which can be accessed only when 
using supercomputer systems. 

An earlier study has shown that the resources of a computing system using only the CPU are not enough to solve 
such scientific problems [3]. Increasing the GPU power and video memory made it possible to use video adapter 
resources for calculations [4]. The GPU utilization depends on the application of parallel algorithms to solve 
computationally intensive problems of aquatic ecology [5—7]. To partially solve the problems of lack of memory and 
computing power on workstations, you can install additional video adapters in PCI-E X16 slots directly and in PCI-E 
X1 slots using PCI-E X1—PCI-E X16 adapters. Thus, the number of video adapters installed on one workstation can be 
increased to 12 [8-11]. 

Heterogeneous computer systems that provide sharing CPU and GPU resources are becoming increasingly popular 
in the scientific community. Application of such systems makes it possible to reduce the computation time of scientific 
problems [12-14]. However, the utilization of a heterogeneous computing environment involves the modernization of 
mathematical models, algorithms and programs that implement them numerically. A heterogeneous system provides 
organizing the calculation process in parallel mode. At the same time, fundamental differences in the construction of 
software systems using CPU and GPU together should be taken into account. 

Materials and Methods. We describe the proposed mathematical models of the computational system, the 
computational grid, as well as the method of decomposition of the computational domain. 

Let D_ beaset of technical characteristics of a computing system, then 

D=D'UD?UD;, (1) 


where D! — a subset of the characteristics of the central processing units (CPU) of a computing system; D? — a 

subset of the characteristics of video adapters (GPU) of a computing system; D? — a subset of RAM characteristics. 
D! =(d'',d2,d'3,d'4), (2) 

where d':! — total number of CPU cores; d'2 — number of streams simultaneously processed by one CPU core; 


d'3 — clock rate, MHz; d! — CPU bus frequency, MHz. 


D= UD? 


Kepu 


={d? | kgey € Kgpy @? e Dz}, (3) 


Kepu €Kepu 
where K py = Toe Newt — multiple video adapter indices; V,,,, — number of computer system video adapters; 


Keepy — video adapter index. Each video adapter is represented as a tuple 


ae ~ (dz, " ee ) : (4) 
where d ioe — amount of video memory of the video adapter with index k,,,,, GB; di; — number of streaming 
multiprocessors. 

D3= (a3, d3) , (5) 


where d3-! — total amount of RAM, GB; d3:2 — clock rate of RAM, MHz. 


Let S —a set of software streams involved in the computing process, then 


S=S'US?, 
S1={1.,N yg}, (6) 
St= U Sz, Sz, = Poo Ns \, 

Kepu €Kgpu kepu 
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where S! —a subset of program streams implementing the calculation process on the CPU; S2 —a subset of CUDA 
streaming blocks implementing the calculation process on GPU streaming multiprocessors; N,, — number of CPU 
program streams involved; Siow — a subset of CUDA streaming blocks that implement the calculation process on 


GPU streaming multiprocessors with index kgpy3 Kgpy = {I,---»Ngpy } — multiple GPU indices; N,,,, — number of 


GPU in the computing system; N,, _— number of CUDA stream blocks involved that implement the calculation 


process on GPU stream multiprocessors with index Kp, . 


Let E bea set of identifiers of program streams. Then, in order to identify program flows in the computing system, 


we assign tuple e of two elements to each element of the set of program flows S: 


VseS deek: e=(n,,n,), (7) 
where n, — index of a computing device in a computing system; n, — index of the CPU program stream or GPU 
streaming unit. 

_ $0, ses! 
= lee. seS?? (8) 
K,, seéS! 
’ 7 ra 2 . E Sian ; ) 


Let us take the computational domain with the following parameters: /, — characteristic size on axis Ox; /, — on 
axis Oy ; |. — on axis Oz. We compare a uniform computational grid of the following type to the specified area: 


W =x, =th,,y, = jh, 2, = kh, 


i=0,n, -1,j =0,n, -1,k =0,n, -1; (10) 
(n, —1)h, a L,,(n, —1)h, = La, —1)h, ee L}, 
where h,, h,, h, — steps of the computational grid in the corresponding spatial directions; n,, n y> 2, —number of 


nodes of the computational grid in the corresponding spatial directions. 


Then, we represent the set of nodes of the computational grid in the form 


G={g,,.1=0,n,-Lj=0,n,-Lk=0,n, -1}, 


Si, j,k = (X15 j2 24) oe 
where g; ,, — computational grid node. 
The number of nodes of the computational grid N,, is calculated from formula 
Ng =n, -n,°N,. (12) 


By the subsection of the computational grid G4 CG (hereinafter — subsection), we mean a subset of the nodes of 


the computational grid G. 


G= U Gr=f{gh|3 eK, gheG lt, 1 Gi=O, (13) 
ky eKy kek, 
where K, = {Lag a — multiple indices of subsections G* of computational grid G; N, — number of subsections 


G4; K,, ; N,, <N; N —set of natural numbers; k, — index of subsection G* . 
Since G4 CG, then 
Gi ={ gh, .i=0.n, 1.70.0) -Lk=0.n, =i], (14) 
where gi , — hode of subsection k,; sign ~ indicates belonging to the subsection; j — node index of subsection 


k, by coordinate y ; ny — number of nodes in subsection k, by coordinate y. 
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ky 


8ij.k 7 (5 ¥) 24), 


k=l 7 15 
5 = thay =( Eat +7} hye =H, 9) 


where n> — number of nodes by coordinate y of 5, -th subsection. 


Under the block of computational grid G4.* (hereinafter — block), we mean a subset of the computational grid 


nodes of subsection G*. 


Gh= YU Gib ={grn [3k eK, ,, ge eGhk}, Gh =O, (16) 


ky ER by Ky ER by 
where K, ,, = aN a — multiple indices of block G4-% of subsection G4; N,, ,, — number of blocks Gi-% ; 
Ky io Net, CN ky — index of block G4-© of subsection G*. 


Since G4i.% c G4, then 


Giuk = {gis ,i=0,n,-1,j7 =0,n'* -1,k =0,n, =} ; (17) 
where oe — node of block k,,k,; sign * indicates belonging to the block; j — node index of block k,,k, by 


coordinate y ; n° — number of nodes in block k,,k, by coordinate y . 
kyky 
Br = (Ris Io 2)> 


ky -1 Nig ty 


- 18 
x, =ih,, y,=|>X ¥ nb 47 \|-h,, z =hh, o) 
i x Jj jal Byel y y k Zz 


where 1} — number of nodes of block b,,p, . 


By a fragment of the computational grid Gi. (hereinafter — fragment), we mean a subset of the nodes of the 


computational grid of block G«» of subsection G*. 


Give = U Ghisko sks = {ghcbok | dk, ek, ; gihikaks (= Ghitets} 
AyeK,, 
1) Gin =D, 


ky eK ,, 


(19) 


where K, ,, ;, ={1, N } — multiple indices of fragments Ges of block G4-% of subsection Gh; N, ,, ,. — 


Bey ie 
number of fragments Ges 5 Ky yas Nu gg, CN 3 k, — index of fragment Gas of block G4. of subsection G'. 
Each index k, of fragment Gs is assigned a tuple of indices es k, Ms designed to store the fragment coordinates 
in the plane xOz , where k, — fragment index by coordinate x; k, — fragment index by coordinate z . 
k, =k, + K,, +k, (20) 
where k, — fragment index by coordinate x; k, — fragment index by coordinate z . 
Number of fragments Ga.» of block G4-% is calculated from the formula 
K, =K,,-K,,; (21) 
where K,,, — number of fragments along axis Ox; K,, — number of fragments by coordinate z . 


Since Ga-es C Ga , then 


Gone ={8,o,, 7=0,7,-1, J=0,8,-1 &=0,7,-1}, (22) 
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where Bi 55 — fragment node; sign _ indicates belonging to the fragment; ike = fragment node indices by 


coordinates x, z; m,, m, — number of nodes of the computational grid in the fragment by coordinates x, z; 


x? Zz 


/,, |, — fragment dimensions by coordinates x, z. 


834 (hs Yj Z,)s 


ky-l = “ ks 1 ih (23) 
x, -( x i, +7 y,=Jh,, % -( > A, +h 
b=1 b=1 
where n, — number of nodes of b -th fragment. 
We introduce a set of comparisons of the computational grid blocks to program flows M! 
M=U U Mis.) (24) 
Ky eKy \ ky eKy, ~ 
where M7; ,, — element of the set M'. 
Let Mj ,, — mapping block G«» to program stream 5, ,, , then 
Mit, = (Gh Ss, ke ) ? (25) 


where s, ,, €S — program flow, computing block Ga . 


In the process of solving hydrodynamic problems on three-dimensional computational grids of large dimension, 
high-performance computing systems and huge amounts of memory for data storage are needed. The resources of one 
computing device are not enough for computing and storing a three-dimensional computational grid with all its data. To 
solve this problem, various methods of decomposition of computational grids followed by the use of parallel calculation 
algorithms in heterogeneous computing environments are proposed [15]. 

For the decomposition of the computational grid, it is required to take into account the performance of computing 
devices involved in calculations. By performance, we mean the number of nodes of the computational grid calculated 
using a given algorithm per unit of time. 


Assume that all computing devices are used for calculations. Then, the total performance of the computing system 


P, is calculated from the formula 


B=P 


Nepu i ; 
3 a Nat 2. Fee Ne ’ (26) 
where P.,, — performance of a single CPU stream; N,, — number of program streams implementing the calculation 
process on the CPU; , ~~ GPU performance with index b ona single streaming multiprocessor; Nf, — number of 


CUDA streaming blocks implementing the calculation process on GPU streaming multiprocessors. 


Then, the number of nodes of the computational grid n’ in the subsection by coordinate y for each GPU with 


index b can be calculated from the formula 
P® 
nb = =| N,. (27) 


In the process of calculating by formula (27), we get the remainder — a certain number of nodes of the 


computational grid. These nodes will be located in RAM. The number of remaining nodes n’ by coordinate y is 
calculated from the formula: 

n’, =n,- > nb, (28) 

To calculate the number of nodes by coordinate y in the blocks of the computational grid processed by GPU 


streaming multiprocessors, we use the formulas: 
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nop =| —— |, b=1,N?’ -1 
yGT ’ 2 92 2 
N [ -1 (29) 
Nn’ -1 
Wer = np ~ > Moors 
where n’,, — number of nodes by coordinate y in computational grid blocks processed by GPU streaming 


multiprocessors with index b , except for the last block; ’,,,, — number of nodes by coordinate y in the last block of 


the computational grid processed by GPU streaming multiprocessors with index 5. 


To calculate the number of nodes of the computational grid by coordinate y in blocks processed by software 


streams implementing the calculation process on the CPU, we use the formulas 


ncPu 
ee et (30) 
Nyon, = NPY —Nyor (N, -1), 

where 1,., — number of nodes of the computational grid by coordinate y , processed by CPU program streams, 


except the last stream; n — number of nodes of the computational grid by coordinate y, processed by CPU 


yCTL 
program streams, in the last stream. 


Calculate the number of the computational grid fragments by coordinate y : 


Nepu 
NJ=N,+ 2 Nt. (31) 


Let the number of fragments N/ and N/ be specified by coordinates x and z, respectively. Then, the number of 


nodes of the computational grid by coordinate x is calculated using the formulas 


Nf -1 (32) 


where nf — number of nodes of the computational grid by coordinate x in all fragments, except the last fragment; 
nf — number of nodes of the computational grid by coordinate x in the last fragment. 


Similarly, the number of nodes of the computational grid is calculated by coordinate z : 


nf = z ; 
* | Nf 1 (33) 
nf =n, —n{ (Nf -1), 


where n/ — number of nodes of the computational grid by coordinate z in all fragments, except the last node; n/ — 
number of nodes of the computational grid by coordinate z in the last fragment. 

Let us describe a model of the parallel-pipeline method. Suppose it is necessary to organize a parallel process of 
computing some function F' on M', and the calculations in each fragment G«-s depend on the values in neighboring 
fragments, each of which has at least one of the indices by coordinates x, y , z, and one less than the current one 
(Fig. 1). 

To organize the parallel-pipeline method, we introduce a set of tuples A, that specify correspondences a between 
the identifiers of program streams e, processing fragments G4..s , to the step numbers of the parallel-pipeline 


method r 


Veek dacAd: a=(e, Ghee), (34) 


where r =1, N, — step number of the parallel-pipeline method; 
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N, — number of steps of the parallel-pipeline method, calculated from the formula 
N = NINT NT 1, (35) 
The full load of all calculators in the proposed parallel-pipeline method starts with step Noosrizr = Nf and ends at 
Step Toosrop = NI N/ . At the same time, the total number of steps with a full load of calculators N,p,, will be 
N par = “oostop ~Moosrarr t1= NE NS -NJ +1. (36) 


The calculation time of some function F’ by the parallel-pipeline method is written as 


= 5 max(T,), (37) 


where T, — vector of values of time spent on processing fragments in parallel mode. 


Ny-1 


(0,0,0) |(1,0,.0) | (2,0,0) | G,0,0) (N,-1,0,0) 
eo —p,» 0 r=0 r=1 r=2 r33 f r= Ny-1 


(0.1.0) (1,1,0) | (2,1,0) | (3,1,0) ee) 
ca 1 r=2 r=3 


0,2,0 1,2,0) | (2,2,0 4.2.0 (N,-1,2,0) 
(0,3,0) (1,3,0) (2,3,0) (3,3,0) o8 S 
e3 —p» 3 r=3 as if 


Fig. 1. Parallel-pipelined computing process 


Research Results. The computational experiments were carried out on K-60 high-performance computing system 
of the Keldysh Applied Mathematics Institute, RAS. A GPU section was used, each node of which was equipped with 
two Intel Xeon Gold 6142 v4 processors, four Nvidia Volta GV100GL video adapters and 768 GB of RAM. 

The computational experiment consisted of two stages — preparatory and basic. At the preparatory stage, the 
correctness of the decomposition of the computational domain into subsections, blocks and fragments was checked by 
step-by-step comparison of values in the nodes of the initial grid and in fragments obtained as a result of decomposition. 
Then, the operation of the flow control algorithm during which the time spent on calculating 1, 8, 16 and 32 fragments 


of the computational grid with a dimension of 50 nodes by spatial coordinates x, y, z, and the same number of CPU 
streams N,, was checked by the iterative alternating-triangular method in parallel mode. Ten repetitions were 
performed with the calculation of the arithmetic mean T, and standard deviation o. Based on the data obtained, 
time T!=T,/N,, spent by each stream on processing one fragment of the computational grid and acceleration 
E=T'(N,)/T!(), equal to the ratio of the processing time T!(N,) of one fragment N,, by streams to the 
corresponding processing time by one stream T!(1) was calculated. The experimental data are given in Table 1. The 


experiment showed that the standard deviation had the smallest value in the case of using 32 parallel CPU streams and 
was 0.026 ms, i.e., using 32 parallel CPU streams when calculating 32 fragments of the computational grid gave a more 
uniform time load of program streams, which generally increased the efficiency of the computing node. At the same 
time, the average value of the calculation of one fragment was 4.14 ms. The dependence of acceleration F on the 


number of streams turned out to be linear E = 0.603+0.804N,, , with a coefficient of determination equal to 0.99. We 
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have found that with an increase in the number of streams, the acceleration of the developed algorithm increases. This 


indicates the efficient use of the subsystem when working with memory. 


Table 1 
Results of the preparatory stage of the computational experiment 


Ng max(T,) , ms o,ms T,, ms E 
1 3.38 0.141 3.38 1.00 
8 3.66 0.042 0.46 7.39 
16 3.94 0.028 0.25 13.73 
32 4.14 0.026 0.13 26.13 


At the basic stage of the computational experiment, a three-dimensional computational domain having dimensions 
of 1,600; 1,600; 200 by spatial coordinates x, y and z, accordingly, was divided into 32 fragments of 50 nodes for 
each of the coordinates x and y. The division into fragments by coordinate z is given in Table 2. For each 
decomposition option with a tenfold repetition, the processing time of the entire computational grid was measured by 


the proposed parallel-conveyor method, and its average value T 


ym Was calculated. Acceleration E,,,, was calculated as 


ratio T,,, to time T,,, of the calculation by sequential version of the algorithm, equal to 6.963 ms. Regression equation 
Ey = 7.3541.97In(N/) with a determination coefficient equal to 0.94 was obtained. Analysis of the results of the 
basic stage of the computational experiment showed a significant slowdown in growth E,,, at N/ >10. Therefore, we 


conclude that splitting into fragments by coordinate z by an amount not exceeding 10 is optimal. 


Table 2 
Results of the main stage of the computational experiment 
Nf nf T,, > Ms Em 
1 200 1033.20 6.74 
2 100 779.00 8.94 
4 50 651.90 10.68 
8 25 588.35 11.84 
20 10 550.22 12.66 


Discussion and Conclusion. As a result of the conducted research, a model of a parallel-pipeline computing process 
was developed by the example of one of the most intensive stages of solving a system of grid equations by a modified 
alternating-triangular iterative method. Its construction was based on decomposition models of a three-dimensional 
uniform computational grid, taking into account the technical characteristics of the equipment used in the calculations. 

The results obtained under the computational experiments validated the effectiveness of the developed method. The 
correctness of the decomposition of the computational domain into subsections, blocks and fragments was also 
confirmed. The operation of the flow control algorithm was verified. At the same time, it was revealed that the standard 
deviation had the smallest value in the case of using 32 parallel CPU streams and is 0.026 ms, i.e., using 32 parallel 
CPU streams when calculating 32 fragments of the computational grid gave a more uniform time load of program 
streams. Here, the average value of the calculation of one fragment was 4.14 ms. 

The results of processing the measurements of the calculation time by the proposed parallel-conveyor method 
showed a significant slowdown in the growth of acceleration when divided into fragments by coordinate z at N/ >10. 
It was found that splitting into fragments by coordinate z by an amount not exceeding 10 was optimal. 
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